RRST : Analysis - sample excluding outlier

Author

Natascha Stoffel

Published

May 23, 2025

0) Data and Preparation

Respiratory Interoception; Analysis of participants of site CH from 15.10.23 until 10.09.2024;

# remove everything that we stored in the environment
rm(list=ls()) 

# Load necessary libraries
library(readr)
library(dplyr)
library(rstatix)
library(gtsummary)
library(ggplot2)
library(ggpubr)
library(openxlsx)
library(gt)
library(stats)
library(grDevices)
library(car)
library(stringr)  
library(lubridate)
library(Hmisc)
library(writexl)
library(tidyr)

### 0.0 PREPARATION -----------------------------------------------------------
# load prepared and cleaned file from 00_preparation: df_yyyymmdd.csv
setwd("/Users/nataschastoffel/Documents/GitHub/interoception_NS")
df_RRST <- read.csv("/Users/nataschastoffel/Documents/GitHub/interoception_NS/data/processed/05-2025/Data.RRST_FINAL_20250509.csv")

# factor variables
factor_variables <- c('group', 'sex', 'psychotropic_medication')
# Apply as.numeric to the selected columns
df_RRST[factor_variables] <- lapply(df_RRST[factor_variables], as.factor)

# numeric variables
numeric_variables <- c('age', 'bmi', 'bdi', 'stai_s', 'stai_t', 'maia_total',
                     'sdq', 'ctq_total', 'ias', 'rrst_sensitivity', 'rrst_slope', 'average_RT', 'metascore',
                     'unpleasantness', 'asthma', 'dizziness', 'breathlessness',
                     'correct_RT', 'incorrect_RT')
# Apply as.numeric to the selected columns
df_RRST[numeric_variables] <- lapply(df_RRST[numeric_variables], as.numeric)

df <- df_RRST
df$maia_commonfactor <- rowSums(df[, c("maia_note", "maia_attreg", "maia_aware", 
                                       "maia_sfreg", "maia_body", "maia_trust")])

1) Having a Look at the Data

1.1) Identifying the Outlier to exclude

# Testing Outlier using boxplot (standard deviation)
df %>% 
  group_by(group) %>%
  identify_outliers(rrst_sensitivity)
# A tibble: 1 × 85
  group pcode date         age sex   handedness diagnosis medication.x smoke_crf
  <fct> <int> <chr>      <dbl> <fct>      <int>     <int>        <int>     <int>
1 1        63 08.05.2024    38 1              1         5            1         0
# ℹ 76 more variables: drugs_crf <int>, bmi <dbl>, bdi <dbl>, stai_t <dbl>,
#   stai_s <dbl>, ias <dbl>, maia_note <dbl>, maia_distr <dbl>,
#   maia_worry <dbl>, maia_attreg <dbl>, maia_aware <dbl>, maia_sfreg <dbl>,
#   maia_body <dbl>, maia_trust <dbl>, maia_total <dbl>, sdq <dbl>,
#   contraception <int>, menopause <int>, mc_phase <int>, sfmdrs <int>,
#   seiz <int>, cgi <int>, ctq_emoab <int>, ctq_physab <int>, ctq_sexab <int>,
#   ctq_emoneg <int>, ctq_physneg <int>, ctq_minim <int>, ctq_total <dbl>, …
# one outlier identified = P063 for rrst_sensitivity 

#### Box-Plot with Jitter
plotBoxplotGroups <- function(df, var , tit ){ 
  plot <- ggplot(df, aes(x = group, y = var, fill = group)) + 
    geom_boxplot(show.legend = TRUE) + 
    labs(y = "Score") + 
    theme_classic() + 
    scale_fill_manual(values = c("#868686FF", "#BB4038")) + 
    scale_x_discrete(labels = c("HC", "FND")) +  
    guides(fill = "none") + 
    ggtitle(tit) +
    xlab("") +
    geom_jitter(shape = 16, position = position_jitter(0.2))
  
  return(plot)
}

p1 <- plotBoxplotGroups(df, (df$rrst_sensitivity), "Respiration Sensitivity") + 
  stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p1

# Filter the rows where rrst_sensitivity is greater than 8 and display pcode (visually detected as extreme)
filtered_pcodes <- df[df$rrst_sensitivity > 8, "pcode"]
print(filtered_pcodes) 
[1] 63
# P063

# Calculate mean and standard deviation of rrst_sensitivity
mean_rrst <- mean(df$rrst_sensitivity)
sd_rrst <- sd(df$rrst_sensitivity)
# Define the threshold for filtering (mean + 2.5 * SD)
threshold <- mean_rrst + 2.5 * sd_rrst # 7.07
# Filter the pcodes for rrst_sensitivity greater than the threshold
filtered_pcodes <- df[df$rrst_sensitivity > threshold, "pcode"]
# View the result
filtered_pcodes
[1] 63
# P063

df_exclusion <- df[df$pcode != "63", ]

p1b <- plotBoxplotGroups(df_exclusion, (df_exclusion$rrst_sensitivity), "Respiration Sensitivity") + 
  stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
print(p1b)

p1b

# Summary of Change - depending on Exclusion of Outlier
summary(df$rrst_sensitivity) 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.458   3.280   3.282   4.350   8.752 
summary(df_exclusion$rrst_sensitivity) 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.429   3.277   3.222   4.307   6.242 
# for statistical reasons we EXCLUDE P063 from the dataset
df <- df[df$pcode != "63", ]

1.2) Testing Distribution of Variables

# Normality of Data -------------------
hist(df$rrst_sensitivity) 
shapiro.test(df$rrst_sensitivity) #  normally distributed

hist(df$rrst_slope) # extremely right skewed
df$rrst_slope_log <- log(df$rrst_slope + 1)  # Add 1 to avoid log(0) and log to normalize skewed variables
hist(df$rrst_slope_log) # still skewed but better (for visualization)
shapiro.test(df$rrst_slope_log) #  NOT normally distributed
shapiro.test(df$rrst_slope_log[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$rrst_slope_log[df$group == "0"]) # controls within themselves are normally distributed

hist(df$metascore) 
shapiro.test(df$metascore)# NOT normally distributed
shapiro.test(df$metascore[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$metascore[df$group == "0"]) # controls within themselves are normally distributed

df$average_RT <- as.numeric(df$average_RT)
hist(df$average_RT) 
shapiro.test(df$average_RT) #  NOT normally distributed

### other variables
shapiro.test(df$age) # normally distributed

shapiro.test(df$bdi) #  NOT normally distributed : p < 0.05
shapiro.test(df$bdi[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$bdi[df$group == "0"]) # controls within themselves are NOT normally distributed

shapiro.test(df$stai_t) #  normally distributed 

shapiro.test(df$stai_s) #  NOT normally distributed 
shapiro.test(df$stai_s[df$group == "1"]) # patients within themselves are normally distributed
shapiro.test(df$stai_s[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05

shapiro.test(df$ctq_total) #  NOT normally distributed : p < 0.05
shapiro.test(df$ctq_total[df$group == "1"]) # patients within themselves are also NOT normally distributed
shapiro.test(df$ctq_total[df$group == "0"]) # controls within themselves are also NOT normally distributed

shapiro.test(df$sdq) #  NOT normally distributed : p < 0.05
shapiro.test(df$sdq[df$group == "1"]) # patients within themselves are NOT normally distributed
shapiro.test(df$sdq[df$group == "0"]) # controls within themselves are NOT normally distributed : p < 0.05

shapiro.test(df$maia_total) # normally distributed 
shapiro.test(df$maia_commonfactor) # normally distributed 
shapiro.test(df$maia_note) # NOT normally distributed 

shapiro.test(df$ias) # NOT normally distributed 

shapiro.test(df$bmi) # NOT normally distributed with p < 0.05
shapiro.test(df$bmi[df$group == "1"]) # within patients bmi is normally distributed 
shapiro.test(df$bmi[df$group == "0"]) # within control bmi is NOT normally distributed with p < 0.05

1.3) Testing Collinearity of Variables

## testing co-dependency of variables by creating a correlation matrix -----
corr_matrix <- cor(df[, c("age","bmi", "bdi", "stai_s", "stai_t", "ctq_total",
                          "ias", "maia_total", "sdq", "rrst_sensitivity", 
                          "rrst_slope", "metascore")], use = "complete.obs")

# Set the cutoff for high correlations, that we then want to exclude
cutoff <- 0.7

# Find the indices of correlations that are greater than 0.7 and exclude self-correlations (diagonal)
high_corr <- which(abs(corr_matrix) > cutoff & abs(corr_matrix) < 1, arr.ind = TRUE)

# Display the pairs of variables with high correlations
high_corr_pairs <- data.frame(
  Var1 = rownames(corr_matrix)[high_corr[, 1]],
  Var2 = colnames(corr_matrix)[high_corr[, 2]],
  Correlation = corr_matrix[high_corr]
)

# View the pairs
print(high_corr_pairs) # stai and bdi are highly correlated, which we do accept though as this is a common clinical happening
cor.test(df$bdi, df$stai_t) #

df <- df %>%
  mutate(anx_dep_SUM = stai_t + bdi) # create one variable for the trait anxiety and depression to use as one variable in the model

2) Group Differences

### Group Differences as Baseline of Variables of Interst
# parametric t-test for NORMALLY distributed data
res<-t.test(df$rrst_sensitivity[df$group=="1"], df$rrst_sensitivity[df$group=="0"])
res # SIG difference between group

    Welch Two Sample t-test

data:  df$rrst_sensitivity[df$group == "1"] and df$rrst_sensitivity[df$group == "0"]
t = -2.1801, df = 73.97, p-value = 0.03243
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.2343241 -0.0554876
sample estimates:
mean of x mean of y 
 2.881984  3.526890 
res<-t.test(df$stai_t[df$group=="1"], df$stai_t[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Welch Two Sample t-test

data:  df$stai_t[df$group == "1"] and df$stai_t[df$group == "0"]
t = 3.9112, df = 82.97, p-value = 0.0001871
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  4.563215 14.006552
sample estimates:
mean of x mean of y 
 46.53488  37.25000 
res<-t.test(df$maia_total[df$group=="1"], df$maia_total[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Welch Two Sample t-test

data:  df$maia_total[df$group == "1"] and df$maia_total[df$group == "0"]
t = -3.679, df = 75.258, p-value = 0.0004374
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -6.547208 -1.947651
sample estimates:
mean of x mean of y 
 19.87420  24.12163 
# non paramatric wilcox test for NOT normally distributed data 
res<-wilcox.test(df$metascore[df$group=="1"], df$metascore[df$group=="0"])
res # no sig difference between groups

    Wilcoxon rank sum exact test

data:  df$metascore[df$group == "1"] and df$metascore[df$group == "0"]
W = 1015, p-value = 0.8962
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$rrst_slope[df$group=="1"], df$rrst_slope[df$group=="0"])
res # no sig difference between groups (p = 0.0629)

    Wilcoxon rank sum exact test

data:  df$rrst_slope[df$group == "1"] and df$rrst_slope[df$group == "0"]
W = 1267, p-value = 0.06209
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$rrst_slope_log[df$group=="1"], df$rrst_slope_log[df$group=="0"])
res # still same result; confirming the data is the same 

    Wilcoxon rank sum exact test

data:  df$rrst_slope_log[df$group == "1"] and df$rrst_slope_log[df$group == "0"]
W = 1267, p-value = 0.06209
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Wilcoxon rank sum test with continuity correction

data:  df$ias[df$group == "1"] and df$ias[df$group == "0"]
W = 734.5, p-value = 0.01813
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$maia_note[df$group=="1"], df$maia_note[df$group=="0"])
res # no sig difference between groups

    Wilcoxon rank sum test with continuity correction

data:  df$maia_note[df$group == "1"] and df$maia_note[df$group == "0"]
W = 986.5, p-value = 0.7191
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$bdi[df$group=="1"], df$bdi[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Wilcoxon rank sum test with continuity correction

data:  df$bdi[df$group == "1"] and df$bdi[df$group == "0"]
W = 1770, p-value = 4.216e-09
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$stai_s[df$group=="1"], df$stai_s[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Wilcoxon rank sum test with continuity correction

data:  df$stai_s[df$group == "1"] and df$stai_s[df$group == "0"]
W = 1585.5, p-value = 1.079e-05
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$ctq_total[df$group=="1"], df$ctq_total[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Wilcoxon rank sum test with continuity correction

data:  df$ctq_total[df$group == "1"] and df$ctq_total[df$group == "0"]
W = 1434.5, p-value = 0.001385
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$sdq[df$group=="1"], df$sdq[df$group=="0"])
res # SIG difference between groups with p < 0.05

    Wilcoxon rank sum test with continuity correction

data:  df$sdq[df$group == "1"] and df$sdq[df$group == "0"]
W = 1838.5, p-value = 1.379e-10
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$ias[df$group=="1"], df$ias[df$group=="0"])
res # no sig difference between groups

    Wilcoxon rank sum test with continuity correction

data:  df$ias[df$group == "1"] and df$ias[df$group == "0"]
W = 734.5, p-value = 0.01813
alternative hypothesis: true location shift is not equal to 0
res<-wilcox.test(df$bmi[df$group=="1"], df$bmi[df$group=="0"])
res # no sig difference between groups

    Wilcoxon rank sum test with continuity correction

data:  df$bmi[df$group == "1"] and df$bmi[df$group == "0"]
W = 1326.5, p-value = 0.01943
alternative hypothesis: true location shift is not equal to 0
#### get mean/median for the main outcome variables per group
mean(df$rrst_sensitivity[df$group=="1"])
[1] 2.881984
mean(df$rrst_sensitivity[df$group=="0"])
[1] 3.52689
sd(df$rrst_sensitivity[df$group=="1"])
[1] 1.619366
sd(df$rrst_sensitivity[df$group=="0"])
[1] 1.128228
median(df$rrst_slope[df$group=="1"])
[1] 63.2525
median(df$rrst_slope[df$group=="0"])
[1] 32.4818
IQR(df$rrst_slope[df$group=="1"])
[1] 2880.422
IQR(df$rrst_slope[df$group=="0"])
[1] 49.09485
median(df$rrst_slope_log[df$group=="1"])
[1] 4.162821
median(df$rrst_slope_log[df$group=="0"])
[1] 3.510839
IQR(df$rrst_slope_log[df$group=="1"])
[1] 4.884927
IQR(df$rrst_slope_log[df$group=="0"])
[1] 1.484032
median(df$metascore[df$group=="1"])
[1] 0.5852489
median(df$metascore[df$group=="0"])
[1] 0.5727064
IQR(df$metascore[df$group=="1"])
[1] 0.07714244
IQR(df$metascore[df$group=="0"])
[1] 0.05455083

2.1) Demographic Summary Table

# define required subgroups
df_female <- df[df$sex == "1", ]
df_male <- df[df$sex == "2", ]

df_HC <- df[df$group == "0", ]
df_FND <- df[df$group == "1", ]


# overview across group and sex distribution
tbl <- df %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male"))) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))
tbl1<-table(tbl$group,tbl$sex)
tbl1
     
      female male
  HC      35   13
  FND     31   12
# overview for mean age across group
summary_stats <- tbl %>%
  group_by(group) %>%
  get_summary_stats(age, type = "full") 
print(summary_stats)
# A tibble: 2 × 14
  group variable     n   min   max median    q1    q3   iqr   mad  mean    sd
  <fct> <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 HC    age         48    18    67     37  28      50  22    14.1  37.8  13.0
2 FND   age         43    18    64     39  31.5    47  15.5  11.9  38.5  11.8
# ℹ 2 more variables: se <dbl>, ci <dbl>
res<-t.test(df$age[df$group == "0"], df$age[df$group == "1"]) # 
res # no sig difference between the groups in age

    Welch Two Sample t-test

data:  df$age[df$group == "0"] and df$age[df$group == "1"]
t = -0.25959, df = 88.975, p-value = 0.7958
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.849281  4.497537
sample estimates:
mean of x mean of y 
 37.81250  38.48837 
# SUMMARY TABLE DEMOGRAPHICS -----
df_tbl <- df %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("male", "female"))) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND"))) %>%
  mutate(psychotropic_medication = factor(psychotropic_medication, levels = c("0", "1"), labels = c("no", "yes"))) 

# Select variable sfor table
Data.t <- select(df_tbl, c("group", "sex", "age", "psychotropic_medication", "smoke_crf", "bmi", "bdi", "stai_t", "ctq_total", "sdq"))

# Create the summary table
tbl1 <- Data.t %>%
    tbl_summary(
    by = group,
    missing = "no",
    statistic = list(
      age = "{mean} ({sd})",
      sex = "{n} ({p})",
      psychotropic_medication = "{n} ({p})",
      smoke_crf = "{n} ({p})",
      bmi = "{median} ({p25}, {p75})",
      bdi = "{median} ({p25}, {p75})",
      stai_t = "{mean} ({sd})",
      ctq_total = "{median} ({p25}, {p75})",
      sdq = "{median} ({p25}, {p75})"
    ),
    digits = list(age = 1, sex = 1, psychotropic_medication = 1, smoke = 1, bmi = 1, bdi = 1, stai_t=1, ctq_total=1, sdq=1),
    label = list(age = "Age, mean (SD)",
                 sex = "Sex, count (%)",
                 psychotropic_medication= "Intake of psychotropic Medication, count (%)",
                 smoke_crf = "Smoking, count (%)",
                 bmi = "Body Mass Index kg/m², median (IQR)",
                 bdi = "Depression using BDI-II, median (IQR), ",
                 stai_t = "Trait Anxiety using STAI, median (IQR)",
                 ctq_total = "Childhood Trauma, using CTQ total, median (IQR)",
                 sdq = "Dissociation using SDQ-20, median (IQR)")
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Demographics Overview**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_p( test = list(
        age ~ "t.test",            # t test for normally distributed variables
        bmi ~ "wilcox.test",      # wilcos test for not-normally distributed variables     
        sex ~ "chisq.test",        # chi-square test for factor variables
        psychotropic_medication ~ "chisq.test",  # Fisher's exact test for small samples
        smoke_crf ~ "chisq.test",  # Fisher's exact test for small samples
        bdi ~ "wilcox.test",
        stai_t ~ "t.test",
        ctq_total = "wilcox.test",
        sdq="wilcox.test"
      )) %>% 
  add_overall(last = FALSE)

# Display the table
tbl1
Demographics Overview
Variable Overall, N = 911 Group p-value2
HC, N = 481 FND, N = 431
Sex, count (%)


>0.9
    male 66.0 (72.5) 35.0 (72.9) 31.0 (72.1)
    female 25.0 (27.5) 13.0 (27.1) 12.0 (27.9)
Age, mean (SD) 38.1 (12.4) 37.8 (13.0) 38.5 (11.8) 0.8
Intake of psychotropic Medication, count (%) 21.0 (23.6) 2.0 (4.3) 19.0 (44.2) <0.001
Smoking, count (%) 15 (16) 6 (13) 9 (21) 0.4
Body Mass Index kg/m², median (IQR) 23.7 (21.7, 26.3) 22.7 (21.7, 25.1) 25.3 (21.8, 27.8) 0.019
Depression using BDI-II, median (IQR), 8.0 (3.0, 15.0) 4.0 (1.0, 8.0) 15.0 (9.5, 23.5) <0.001
Trait Anxiety using STAI, median (IQR) 41.6 (12.1) 37.3 (10.3) 46.5 (12.1) <0.001
Childhood Trauma, using CTQ total, median (IQR) 37.0 (30.0, 54.5) 34.0 (29.0, 41.3) 46.0 (35.0, 60.5) 0.001
Dissociation using SDQ-20, median (IQR) 27.0 (22.0, 36.0) 23.0 (21.0, 26.0) 36.0 (28.5, 45.0) <0.001
1 n (%); Mean (SD); Median (IQR)
2 Pearson’s Chi-squared test; Welch Two Sample t-test; Wilcoxon rank sum test
    # save table as excel
    tbl1_df <- tbl1 %>% as_tibble()
    write_xlsx(tbl1_df, "2.1_Table_demographics.xlsx")

2.2) FND Summary Table

# SUMMARY TABLE for FND symptoms -----
# create variable of duation of months

df_FND <- df_FND %>%
  mutate(duration_months = interval(date_symptom_onset, date) %/% months(1))
as.numeric(df_FND$duration_months)
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
[26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
# select variables for table
Data.t <- dplyr::select(df_FND, c("sex", "duration_months", "fds", "motor", "weakness", "sensory", "pppd", "cognitive"))

Data.t <- Data.t  %>%
  mutate(sex = factor(sex, levels = c("1", "2"), labels = c("female", "male")))

# Convert 0/1 values to "no"/"yes" for table, so that it will only show the count "yes"
Data.t <- Data.t %>%
  mutate(
    across(c(fds, motor, weakness, sensory, pppd, cognitive), 
           ~ factor(.x, levels = c(0, 1), labels = c("no", "yes")))
  ) %>%
  dplyr::select(sex, duration_months, fds, motor, weakness, sensory, pppd, cognitive)

# Create the summary table (which symptoms present with how many patients ; per patient more than one symptom-category possible))
tbl1 <- Data.t %>%
  tbl_summary(
    by = sex,
    missing = "no",
    statistic = list(
      duration_months = "{mean} ({sd})",
      fds = "{n} ({p})",
      motor = "{n} ({p})",
      weakness = "{n} ({p})",
      sensory = "{n} ({p})",
      pppd = "{n} ({p})",
      cognitive = "{n} ({p})"
    ),
    digits = list(medication = 1, duration_months = 1, fds = 1, motor = 1, weakness = 1, sensory = 1, pppd = 1, cognitive = 1),
    label = list(
      sex = "Sex",
      duration_months = "Symptom Duration [in months]",
      fds = "Functional Dissociative Seizures [yes]",
      motor = "motor + symptoms [yes]",
      weakness = "motor - symptoms [yes]",
      sensory = "sensory symptoms [yes]",
      pppd = "dizziness (PPPD) [yes]",
      cognitive = "cognitive symptoms [yes]"
    )
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Clinical Characteristics**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_p(
    test = list(
      duration_months = "t.test", 
      fds = "fisher.test", 
      motor = "fisher.test",
      weakness = "fisher.test", 
      sensory = "fisher.test", 
      ctq_total = "fisher.test", 
      pppd = "fisher.test", 
      cognitive = "fisher.test"
    )  
    ) %>%   add_overall(last = FALSE)

tbl1
Clinical Characteristics
Variable Overall, N = 431 Group p-value2
female, N = 311 male, N = 121
Symptom Duration [in months] NA (NA) NA (NA) NA (NA)
Functional Dissociative Seizures [yes] 8.0 (18.6) 8.0 (25.8) 0.0 (0.0) 0.082
motor + symptoms [yes] 20.0 (46.5) 13.0 (41.9) 7.0 (58.3) 0.5
motor - symptoms [yes] 29.0 (67.4) 20.0 (64.5) 9.0 (75.0) 0.7
sensory symptoms [yes] 19.0 (44.2) 15.0 (48.4) 4.0 (33.3) 0.5
dizziness (PPPD) [yes] 2.0 (4.7) 2.0 (6.5) 0.0 (0.0) >0.9
cognitive symptoms [yes] 4.0 (9.3) 4.0 (12.9) 0.0 (0.0) 0.6
1 Mean (SD); n (%)
2 Fisher’s exact test
  # save table as excel
  tbl1_df <- tbl1 %>% as_tibble()
  write_xlsx(tbl1_df, "2.2_Table_clinical_characteristics.xlsx")

2.3) Visualization of Group Differences

Here we visualize the different distribution per group across all variables

df$group <- as.factor(df$group)

# plot function for  normally dsitributed data, allowing the use of the t test (and the mean)
plotViolinGroups_mean <- function(df, var , tit ){ 
  plot <-ggplot(df, aes(x=group, y=var, fill=group)) + 
    geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8)+  
    stat_summary(show.legend = FALSE, fun = mean, geom = "point", color = "black", size = 2) +  
    labs(y = "Score") + 
    theme_classic()+ 
    scale_fill_manual(values = c("#868686FF", "#BB4038"))+ 
    scale_x_discrete(labels=c("Controls", "Patients"))+  
    guides(fill=guide_legend(title="Timepoint"))+ 
    ggtitle(tit) +
    xlab("")
  
  return(plot)
}

# plot function for non normally dsitributed data; the wilcox test (and the median)
plotViolinGroups_median <- function(df, var , tit ){ 
  plot <-ggplot(df, aes(x=group, y=var, fill=group)) +  
    geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8)+ 
    stat_summary(show.legend = FALSE, fun = median, geom = "point", color = "black", size = 2) + 
    labs(y = "Score") + 
    theme_classic()+ 
    scale_fill_manual(values = c("#868686FF", "#BB4038"))+ 
    scale_x_discrete(labels=c("Controls", "Patients"))+  
    ggtitle(tit) +
    xlab("")
  
  return(plot)
}

# adding jitter points for visualization in publication
jitteredViolin <- function(df, var, tit){  
  ggplot(df, aes(x = group, y = {{ var }}, fill = group)) +  
    geom_violin(show.legend = FALSE, width = 0.7, trim = FALSE, alpha = 0.8) +  
    geom_jitter(width = 0.1, size = 1.2, color = "black", alpha = 0.6) +  
    stat_summary(show.legend = FALSE, fun = median, geom = "point", color = "black", size = 2) +  
    labs(y = "Score") +  
    theme_classic() +  
    theme(legend.position = "none") +
    scale_fill_manual(values = c("#868686FF", "#BB4038")) +  
    scale_x_discrete(labels = c("Controls", "Patients")) +  
    ggtitle(tit) +  
    xlab("")  
}

p1 <- plotViolinGroups_median(df,df$bdi,"BDI - Depression")+ ggpubr::stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p1

  pdf("2.3_Plot_BDI_violin.pdf")
  p1
  dev.off()
quartz_off_screen 
                2 
p2 <- plotViolinGroups_mean(df,df$stai_s,"STAI - State Anxiety")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p3 <- plotViolinGroups_median(df,df$stai_t,"STAI - Trait Anxiety")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p2, p3, ncol = 2)

  pdf("2.3_Plot_STAI_violin.pdf")
  gridExtra::grid.arrange(p2, p3, ncol = 2)
  dev.off()
quartz_off_screen 
                2 
p4 <- plotViolinGroups_median(df,df$ctq_total,"CTQ - Childhood Trauma")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p4

  pdf("2.3_Plot_CTQ_violin.pdf")
  p4
  dev.off()
quartz_off_screen 
                2 
p5 <- plotViolinGroups_median(df,df$sdq,"SDQ-20 - Dissociation")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p5

    pdf("2.3_Plot_SDQ-20_violin.pdf")
  p5
  dev.off()
quartz_off_screen 
                2 
p6 <- plotViolinGroups_mean(df,df$maia_total,"MAIA - total score")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
  
p6common <- plotViolinGroups_mean(df,df$maia_commonfactor,"MAIA - common factor")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
  
p6noticing <- plotViolinGroups_median(df,df$maia_note,"MAIA subscale - Noticing")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)

  gridExtra::grid.arrange(p6, p6common, p6noticing, ncol = 3)

  pdf("2.3_Plots_MAIA_violin.pdf")
  gridExtra::grid.arrange(p6, p6common, p6noticing, ncol = 3)
  dev.off()
quartz_off_screen 
                2 
p7 <- plotViolinGroups_median(df,df$ias,"IAS - Interoceptive Accuracy")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p7

  pdf("2.3_Plot_IAS_violin.pdf")
  p7
  dev.off()
quartz_off_screen 
                2 
p8 <- plotViolinGroups_mean(df,df$rrst_sensitivity,"Respiratory Sensitivity (RRST alpha)")+ stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
p9 <- plotViolinGroups_median(df,df$rrst_slope_log,"Respiratory Slope (RRST beta)")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
p10 <- plotViolinGroups_median(df,df$metascore,"Respiratory Metacognition")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
gridExtra::grid.arrange(p8, p9, p10, ncol = 2)

  pdf("2.3_Plots_RRST_violin.pdf")
  gridExtra::grid.arrange(p8, p9, p10,  ncol = 2)
  dev.off()
quartz_off_screen 
                2 
  Figure2c <- jitteredViolin(df,df$rrst_sensitivity,"Respiratory Sensitivity (RRST alpha)") +
  stat_compare_means(method = "t.test", paired = FALSE, label.x = 1.3)
  Figure2c

  Figure2d <- jitteredViolin(df,df$rrst_slope_log, "Respiratory Slope (RRST beta)") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
  Figure2d

  Figure3c <- jitteredViolin(df,df$metascore,"Respiratory Metacognition (AUC)") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
  Figure3c

 #### save images as svg to merge into a tiff for publication 
  ggsave("Figure2c.svg", 
       plot = Figure2c, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "svg")
  
ggsave("Figure2d.svg", 
       plot = Figure2d, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "svg")
    
  ggsave("Figure3c.svg", 
       plot = Figure3c, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "svg")

# Table for dependent Variables: group differences
df_tbl <- df %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND"))) 
  
Data.t <- select(df_tbl, c("group", "maia_total", "maia_commonfactor", "maia_note", "ias", "rrst_sensitivity","rrst_slope", "metascore"))
library(gtsummary)
tbl2 <- Data.t %>%
  tbl_summary(
    by = group,
    missing = "no",
    statistic = list(
      maia_total = "{mean} ({sd})",
      maia_commonfactor = "{mean} ({sd})",
      maia_note = "{median} ({IQR})",
      ias = "{median} ({IQR})",
      rrst_sensitivity = "{mean} ({sd})",
      rrst_slope = "{median} ({IQR})",
      metascore = "{median} ({IQR})"
    ),
    digits = list(maia_total = 2, maia_commonfactor = 2, maia_note = 2, ias = 2, 
                  rrst_sensitivity = 2, rrst_slope = 2, metascore = 2),
    label = list(maia_total= "Interoceptive Awareness Self-Report - MAIA, Mean, (SD)",
                  maia_commonfactor= "MAIA common factor, Mean, (SD)",
                 maia_note= "MAIA subscale Noticing, Median, (IQR)",
                 ias= "Interoceptive Accuracy Self-Report - IAS, Median (IQR)",
                 rrst_sensitivity= "Respiratory Sensitivity - RRST, Mean (SD)",
                rrst_slope= "Interoceptive Decision Prescision - RRST slope, Median (IQR)",
                 metascore= "Respiratory Metacognition - AUC, Median (IQR)"
                 )
  ) %>%
  modify_header(list(label ~ "**Variable**")) %>%
  modify_caption("**Interoception Overview**") %>%
  modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") %>%
  add_p(
    test = list(
      maia_total = "t.test", 
      maia_commonfactor = "wilcox.test", 
      maia_note = "t.test",  
      ias = "wilcox.test",        
      rrst_sensitivity = "t.test",
      rrst_slope = "wilcox.test",
      metascore = "wilcox.test")
    ) %>%
    add_overall(last = FALSE)
tbl2
Interoception Overview
Variable Overall, N = 911 Group p-value2
HC, N = 481 FND, N = 431
Interoceptive Awareness Self-Report - MAIA, Mean, (SD) 22.11 (5.78) 24.12 (4.49) 19.87 (6.27) <0.001
MAIA common factor, Mean, (SD) 17.15 (5.49) 19.07 (3.82) 15.01 (6.27) <0.001
MAIA subscale Noticing, Median, (IQR) 3.25 (1.00) 3.25 (0.75) 3.25 (1.63) 0.4
Interoceptive Accuracy Self-Report - IAS, Median (IQR) 85.00 (18.50) 89.00 (14.25) 80.00 (25.50) 0.018
Respiratory Sensitivity - RRST, Mean (SD) 3.22 (1.41) 3.53 (1.13) 2.88 (1.62) 0.032
Interoceptive Decision Prescision - RRST slope, Median (IQR) 39.01 (421.44) 32.48 (49.09) 63.25 (2,880.42) 0.062
Respiratory Metacognition - AUC, Median (IQR) 0.57 (0.06) 0.57 (0.05) 0.59 (0.08) 0.9
1 Mean (SD); Median (IQR)
2 Welch Two Sample t-test; Wilcoxon rank sum test; Wilcoxon rank sum exact test
    tbl2_df <- tbl2 %>% as_tibble()
    write_xlsx(tbl2_df, "2.3_Table_dependent_variables.xlsx")

2.4) Effect Sizes and CI

Here we calculate the effect sizes (cohen’s d or OR) and the confidence intervals of the variables to be reported

# calculation of cohens d for variable of interest

# Define the function
calculate_cohens_d <- function(data, variable, group_var, group1, group2) {
  # Subset the data for each group
  group1_data <- data[[variable]][data[[group_var]] == group1]
  group2_data <- data[[variable]][data[[group_var]] == group2]
  
  # Calculate means and standard deviations
  mean1 <- mean(group1_data, na.rm = TRUE)
  mean2 <- mean(group2_data, na.rm = TRUE)
  sd1 <- sd(group1_data, na.rm = TRUE)
  sd2 <- sd(group2_data, na.rm = TRUE)
  
  # Calculate sample sizes
  n1 <- sum(!is.na(group1_data))
  n2 <- sum(!is.na(group2_data))
  
  # Calculate pooled standard deviation
  sd_pooled <- sqrt(((n1 - 1) * sd1^2 + (n2 - 1) * sd2^2) / (n1 + n2 - 2))
  
  # Calculate Cohen's d
  cohens_d <- (mean1 - mean2) / sd_pooled
  
  # Return Cohen's d
  return(cohens_d)
}

# Calculate Cohen's d for variables of interest
cohens_d_rrst <- calculate_cohens_d(df, "rrst_sensitivity", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_rrst) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC
[1] -0.4666721
cohens_d_slope <- calculate_cohens_d(df, "rrst_slope", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_slope) # cohens d of 0.509 means with an effect size of 0.51 FND has HIGHER values than HC
[1] 0.5088956
cohens_d_meta <- calculate_cohens_d(df, "metascore", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_meta) # cohens d of -0.128 means with an effect size of 0.128 FND has LOWER values than HC
[1] -0.1279049
cohens_d_maia <- calculate_cohens_d(df, "maia_total", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia) # cohens d of -0.786 means with an effect size of 0.786 FND has LOWER values than HC
[1] -0.7864429
cohens_d_maia <- calculate_cohens_d(df, "maia_commonfactor", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia) # cohens d of -0.793 means with an effect size of 0.793 FND has LOWER values than HC
[1] -0.7925178
cohens_d_maia_note <- calculate_cohens_d(df, "maia_note", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_maia_note) # cohens d of -0.182 means with an effect size of 0.18 FND has LOWER values than HC
[1] -0.1824447
cohens_d_ias <- calculate_cohens_d(df, "ias", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_ias) # cohens d of -0.651 means with an effect size of 0.651 FND has LOWER values than HC
[1] -0.6510104
cohens_d_bdi <- calculate_cohens_d(df, "bdi", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_bdi) # cohens d of 1.382615 means with an effect size of 1.382615 FND has HIGHER values than HC
[1] 1.382615
cohens_d_stai_t <- calculate_cohens_d(df, "stai_t", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_stai_t) # cohens d of 0.829 means with an effect size of 0.829 FND has HIGHER values than HC
[1] 0.8285794
cohens_d_stai_s <- calculate_cohens_d(df, "stai_s", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_stai_s) # cohens d of 1.080 means with an effect size of 1.080 FND has HIGHER values than HC
[1] 1.080068
cohens_d_ctq_total <- calculate_cohens_d(df, "ctq_total", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_ctq_total) # cohens d of 0.662 means with an effect size of 0.662 FND has HIGHER values than HC
[1] 0.6621177
cohens_d_sdq <- calculate_cohens_d(df, "sdq", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_sdq) # cohens d of 1.500 means with an effect size of 1.500 FND has HIGHER values than HC
[1] 1.500201
cohens_d_bmi <- calculate_cohens_d(df, "bmi", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_bmi) # cohens d of 0.557 means with an effect size of 0.6621177 FND has HIGHER values than HC
[1] 0.5569856
# for binary variable we calculate the ODDS RATIO
model <- glm(group ~ psychotropic_medication, data = df, family = binomial)
summary(model)

Call:
glm(formula = group ~ psychotropic_medication, family = binomial, 
    data = df)

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -0.6061     0.2538  -2.389 0.016912 *  
psychotropic_medication1   2.8574     0.7854   3.638 0.000275 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 123.28  on 88  degrees of freedom
Residual deviance: 101.51  on 87  degrees of freedom
  (2 observations deleted due to missingness)
AIC: 105.51

Number of Fisher Scoring iterations: 4
# Extract CI and Odds Ratio
confint_logit <- confint(model)  # CI on the log-odds scale
confint_or <- exp(confint_logit) # Exponentiate to get CI for OR
odds_ratio <- exp(coef(model)) # Calculate Odds Ratios
results <- cbind(odds_ratio, confint_or) # Combine OR and CI
results
                         odds_ratio     2.5 %      97.5 %
(Intercept)               0.5454545 0.3268162   0.8883923
psychotropic_medication1 17.4166659 4.5447406 115.3119754

2.5) MAIA subscores

 maia<-dplyr::select(df, pcode, group, maia_attreg, maia_note, maia_distr, maia_worry, maia_aware, maia_sfreg, maia_body, maia_trust)
  
  # Make long data
  maia_long <- gather(maia, key = "subscale", value = "score", -pcode, -group) # make data set long (so that every p code has 6 rows, one for each of the subscores of ctq)
  maia_long$subscale<-as.factor(maia_long$subscale)
 maia_long$subscale <- factor(maia_long$subscale, levels = rev(levels(maia_long$subscale)))
  maia_long$score<-as.numeric(maia_long$score)
  maia_long$group<-as.factor(maia_long$group)
  maia_long$group <- factor(maia_long$group, levels = c(0, 1), labels = c("HC", "FND"))
  
  # Calculate mean, median, CI, SE
  library(dplyr)
  alpha=0.05
  maia.summary <- maia_long %>%
    group_by(group,subscale)%>%
    dplyr::summarise(
      se = sd(score) / sqrt(length(score)),
      t=qt((1-alpha)/2 + .5, length(score)-1),
      CI=t*se, 
      score = mean(score), .groups = 'drop') 

  maia.summary # according to this order name the lables in the following plot p
# A tibble: 16 × 6
   group subscale        se     t    CI score
   <fct> <fct>        <dbl> <dbl> <dbl> <dbl>
 1 HC    maia_worry  0.144   2.01 0.289  2.85
 2 HC    maia_trust  0.161   2.01 0.324  3.51
 3 HC    maia_sfreg  0.151   2.01 0.303  2.82
 4 HC    maia_note   0.0912  2.01 0.183  3.43
 5 HC    maia_distr  0.111   2.01 0.223  2.21
 6 HC    maia_body   0.136   2.01 0.273  2.69
 7 HC    maia_aware  0.119   2.01 0.240  3.68
 8 HC    maia_attreg 0.119   2.01 0.239  2.93
 9 FND   maia_worry  0.166   2.02 0.334  2.97
10 FND   maia_trust  0.224   2.02 0.453  2.39
11 FND   maia_sfreg  0.201   2.02 0.407  1.88
12 FND   maia_note   0.182   2.02 0.367  3.26
13 FND   maia_distr  0.178   2.02 0.359  1.90
14 FND   maia_body   0.207   2.02 0.417  2.10
15 FND   maia_aware  0.214   2.02 0.432  3.08
16 FND   maia_attreg 0.176   2.02 0.355  2.30
  # STATISTICS
  # Multiple t-tests using FDR correction for multiple comparisons
  library(rstatix)
  stat.test <- maia_long %>%
    group_by(subscale) %>%
    t_test(score ~ group) %>%
    adjust_pvalue(method = "fdr") %>% 
    add_significance()
  stat.test # subscales worrying, trusting, body and attention regulation are sig difference across groups
# A tibble: 8 × 11
  subscale    .y.   group1 group2    n1    n2 statistic    df        p    p.adj
  <fct>       <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl>
1 maia_worry  score HC     FND       48    43    -0.555  85.7 0.58     0.58    
2 maia_trust  score HC     FND       48    43     4.08   78.0 0.000109 0.000872
3 maia_sfreg  score HC     FND       48    43     3.73   79.8 0.000353 0.00141 
4 maia_note   score HC     FND       48    43     0.842  62.3 0.403    0.461   
5 maia_distr  score HC     FND       48    43     1.48   71.4 0.144    0.192   
6 maia_body   score HC     FND       48    43     2.40   73.7 0.0189   0.0302  
7 maia_aware  score HC     FND       48    43     2.45   66.4 0.0169   0.0302  
8 maia_attreg score HC     FND       48    43     2.97   75.0 0.00404  0.0108  
# ℹ 1 more variable: p.adj.signif <chr>
  # Filter for significant results (both * and **)
  significant_results <- stat.test %>%
    filter(p.adj.signif %in% c("*", "**")) %>%
    select(subscale, p.adj.signif)  # Keep both subscale and significance level
  
  # VISUALIZATION
  # Plot
  library(ggsignif)
  maia.summary <- maia.summary %>%
    left_join(significant_results, by = "subscale") %>%
    mutate(significance = ifelse(is.na(p.adj.signif), "", p.adj.signif))  # Add a column for asterisks
  
  dodge <- position_dodge(width=0.9)
  maiaplot <- maia.summary %>%
  ggplot(aes(y = score, x = subscale, ymin = score - CI, ymax = score + CI, fill = group)) +
  geom_bar(position = position_dodge(), stat = "identity") +
  geom_errorbar(position = position_dodge(width = 0.9), width = 0.3) +
  labs(x = "MAIA Subscales", y = "Scores") +
  scale_fill_manual(values = c("#868686FF", "#BB4038")) +
  ggtitle("MAIA - Interoceptive Awareness") +
  theme_classic(base_size = 9) +
  scale_x_discrete(labels = c("Not-Worrying", "Trusting", "Self-Regulation", "Noticing", 
                              "Not-Distracting", "Body Listening", "Emotional Awareness", "Attention Regulation")) +
  guides(fill = guide_legend(title = "Groups"))
    maiaplot

    pdf("2.5_MAIAsubscales.pdf")
    maiaplot
    dev.off()
quartz_off_screen 
                2 

3) Testing potential Confounders for RRST sensitivity

3.2) Anxiety

# RRST Sensitivity GROUP
rrst_affective <- lm(formula = rrst_sensitivity ~  group,
                 data = df) 
summary(rrst_affective)  # GROUP

Call:
lm(formula = rrst_sensitivity ~ group, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.88198 -0.94639 -0.07559  0.99151  3.07322 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.5269     0.1995  17.682   <2e-16 ***
group1       -0.6449     0.2902  -2.223   0.0288 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.382 on 89 degrees of freedom
Multiple R-squared:  0.05258,   Adjusted R-squared:  0.04194 
F-statistic:  4.94 on 1 and 89 DF,  p-value: 0.02878
confint(rrst_affective, level = 0.95)  
                2.5 %     97.5 %
(Intercept)  3.130559  3.9232197
group1      -1.221465 -0.0683471
# RRST Sensitivity GROUP controlled for affective symptoms
rrst_affective <- lm(formula = rrst_sensitivity ~  group + anx_dep_SUM,
                 data = df) 
summary(rrst_affective)  # no sig anymore

Call:
lm(formula = rrst_sensitivity ~ group + anx_dep_SUM, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.94169 -0.91149 -0.05828  1.02851  3.07928 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.715201   0.395404   9.396 6.28e-15 ***
group1      -0.552118   0.336279  -1.642    0.104    
anx_dep_SUM -0.004418   0.007999  -0.552    0.582    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.387 on 88 degrees of freedom
Multiple R-squared:  0.05586,   Adjusted R-squared:  0.0344 
F-statistic: 2.603 on 2 and 88 DF,  p-value: 0.07974
confint(rrst_affective, level = 0.95)  
                  2.5 %     97.5 %
(Intercept)  2.92941812 4.50098406
group1      -1.22040198 0.11616628
anx_dep_SUM -0.02031377 0.01147804
# RRST Sensitivity GROUP controlled for affective symptoms
rrst_anxiety <- lm(formula = rrst_sensitivity ~  group + stai_t,
                 data = df) 
summary(rrst_anxiety)  # no sig anymore

Call:
lm(formula = rrst_sensitivity ~ group + stai_t, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.89893 -0.93955 -0.09417  0.97389  3.07027 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.451869   0.529118   6.524 4.21e-09 ***
group1      -0.663605   0.316270  -2.098   0.0388 *  
stai_t       0.002014   0.013144   0.153   0.8786    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.39 on 88 degrees of freedom
Multiple R-squared:  0.05284,   Adjusted R-squared:  0.03131 
F-statistic: 2.454 on 2 and 88 DF,  p-value: 0.09177
confint(rrst_anxiety, level = 0.95)  
                  2.5 %      97.5 %
(Intercept)  2.40035747  4.50337973
group1      -1.29212600 -0.03508498
stai_t      -0.02410788  0.02813585
# RRST Sensitivity affective symptoms
rrst_affective <- lm(formula = rrst_sensitivity ~  anx_dep_SUM,
                 data = df) 
summary(rrst_affective)  # anx_dep sum score by itself NOT sig for rrst_sensitivity

Call:
lm(formula = rrst_sensitivity ~ anx_dep_SUM, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2496 -0.8413  0.0435  1.0227  2.8697 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.799079   0.395808   9.598 2.17e-15 ***
anx_dep_SUM -0.010979   0.006995  -1.570     0.12    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.401 on 89 degrees of freedom
Multiple R-squared:  0.02693,   Adjusted R-squared:  0.016 
F-statistic: 2.463 on 1 and 89 DF,  p-value: 0.1201
confint(rrst_affective, level = 0.95)  
                  2.5 %      97.5 %
(Intercept)  3.01261753 4.585539931
anx_dep_SUM -0.02487721 0.002919802
# RRST Sensitivity anxiety symptoms
rrst_anxiety <- lm(formula = rrst_sensitivity ~  stai_t,
                 data = df) 
summary(rrst_anxiety)  # anx_dep sum score by itself NOT sig for rrst_sensitivity

Call:
lm(formula = rrst_sensitivity ~ stai_t, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2444 -0.7941  0.0841  1.1044  2.9103 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.581435   0.535453   6.689 1.92e-09 ***
stai_t      -0.008629   0.012356  -0.698    0.487    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.416 on 89 degrees of freedom
Multiple R-squared:  0.00545,   Adjusted R-squared:  -0.005725 
F-statistic: 0.4877 on 1 and 89 DF,  p-value: 0.4868
confint(rrst_anxiety, level = 0.95)  
                  2.5 %     97.5 %
(Intercept)  2.51750016 4.64536953
stai_t      -0.03317991 0.01592229

3.2) RRST discomfort

## first we test whether the different control-questions correlate
corr_matrix <- cor(df[, c("breathlessness","dizziness", "unpleasantness", "asthma")], use = "complete.obs")

# Set the cutoff for high correlations, that we then want to exclude
cutoff <- 0.5

# Find the indices of correlations that are greater than 0.7 and exclude self-correlations (diagonal)
high_corr <- which(abs(corr_matrix) > cutoff & abs(corr_matrix) < 1, arr.ind = TRUE)

# Display the pairs of variables with high correlations
high_corr_pairs <- data.frame(
  Var1 = rownames(corr_matrix)[high_corr[, 1]],
  Var2 = colnames(corr_matrix)[high_corr[, 2]],
  Correlation = corr_matrix[high_corr]
)

# View the pairs
print(high_corr_pairs) # no high correlation; no general factor of discomfort
[1] Var1        Var2        Correlation
<0 rows> (or 0-length row.names)
### calculate correlations to respiratory sensitivity (and control for multiple comparision)
# Run all correlations and store p-values
p_values <- c(
  cor.test(df$breathlessness, df$rrst_sensitivity)$p.value,
  cor.test(df$asthma, df$rrst_sensitivity)$p.value,
  cor.test(df$unpleasantness, df$rrst_sensitivity)$p.value,
  cor.test(df$dizziness, df$rrst_sensitivity)$p.value
)

# Apply FDR correction
p_adjusted <- p.adjust(p_values, method = "fdr")

# Optional: associate each result with variable name
results <- data.frame(
  variable = c("breathlessness", "asthma", "unpleasantness", "dizziness"),
  raw_p = p_values,
  p_fdr = p_adjusted)
print(results)
        variable       raw_p       p_fdr
1 breathlessness 0.025679985 0.051359970
2         asthma 0.001684408 0.006737632
3 unpleasantness 0.302098918 0.402798557
4      dizziness 0.921408315 0.921408315
# Calculate Cohen's d for variables of interest
cohens_d_discomfort <- calculate_cohens_d(df, "breathlessness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC
[1] 0.6465991
cohens_d_discomfort <- calculate_cohens_d(df, "asthma", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC
[1] 0.4498372
cohens_d_discomfort <- calculate_cohens_d(df, "unpleasantness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC
[1] 0.4801304
cohens_d_discomfort <- calculate_cohens_d(df, "dizziness", "group", 1, 0) # group 1=FND, group=0 HC
print(cohens_d_discomfort) # cohens d of -0.467 means with an effect size of 0.46 FND has LOWER values than HC
[1] 0.4457795
## testing linear models for the effect on the identified group difference
# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_discomfort <- lm(formula = rrst_sensitivity ~  group + unpleasantness,
                 data = df) 
summary(rrst_discomfort)  # not sig anymore; p = 0.0691

Call:
lm(formula = rrst_sensitivity ~ group + unpleasantness, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.01017 -0.92969 -0.01731  1.06434  2.95843 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     3.665167   0.271946  13.478   <2e-16 ***
group1         -0.558640   0.303399  -1.841   0.0691 .  
unpleasantness -0.003455   0.005862  -0.589   0.5572    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.38 on 85 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.05026,   Adjusted R-squared:  0.02791 
F-statistic: 2.249 on 2 and 85 DF,  p-value: 0.1117
confint(rrst_discomfort, level = 0.95)  
                     2.5 %      97.5 %
(Intercept)     3.12446522 4.205868773
group1         -1.16187795 0.044597149
unpleasantness -0.01510882 0.008199755
# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_asthma <- lm(formula = rrst_sensitivity ~  group + asthma,
                 data = df) 
summary(rrst_asthma)  # group not sig, but asthma with p. = 0.005

Call:
lm(formula = rrst_sensitivity ~ group + asthma, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2401 -0.6751  0.0155  0.9793  2.7620 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.657422   0.195750  18.684  < 2e-16 ***
group1      -0.417360   0.289355  -1.442  0.15287    
asthma      -0.022484   0.007853  -2.863  0.00528 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.32 on 85 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.1303,    Adjusted R-squared:  0.1098 
F-statistic: 6.365 on 2 and 85 DF,  p-value: 0.002656
confint(rrst_asthma, level = 0.95)  
                  2.5 %       97.5 %
(Intercept)  3.26821773  4.046625986
group1      -0.99267486  0.157954761
asthma      -0.03809798 -0.006869975
# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_breath <- lm(formula = rrst_sensitivity ~  group + breathlessness,
                 data = df) 
summary(rrst_breath)  # both factors not sig

Call:
lm(formula = rrst_sensitivity ~ group + breathlessness, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1858 -0.7460 -0.0539  0.9959  2.9155 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     3.701317   0.215022  17.214   <2e-16 ***
group1         -0.437026   0.305481  -1.431   0.1562    
breathlessness -0.010756   0.006226  -1.728   0.0877 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.359 on 85 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.07873,   Adjusted R-squared:  0.05706 
F-statistic: 3.632 on 2 and 85 DF,  p-value: 0.03065
confint(rrst_breath, level = 0.95)  
                     2.5 %      97.5 %
(Intercept)     3.27379618 4.128838668
group1         -1.04440446 0.170352839
breathlessness -0.02313453 0.001621564
# RRST Sensitivity GROUP controlled for discomfort of the task
rrst_dizzy <- lm(formula = rrst_sensitivity ~  group + dizziness,
                 data = df) 
summary(rrst_dizzy)  # group still sig p = 0.042

Call:
lm(formula = rrst_sensitivity ~ group + dizziness, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.99484 -0.79845 -0.04541  0.97905  3.02539 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.507924   0.245147  14.310   <2e-16 ***
group1      -0.624278   0.302629  -2.063   0.0422 *  
dizziness    0.001906   0.005379   0.354   0.7239    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.382 on 85 degrees of freedom
  (3 observations deleted due to missingness)
Multiple R-squared:  0.04778,   Adjusted R-squared:  0.02538 
F-statistic: 2.133 on 2 and 85 DF,  p-value: 0.1248
confint(rrst_dizzy, level = 0.95)  
                   2.5 %      97.5 %
(Intercept)  3.020507382  3.99534123
group1      -1.225986511 -0.02257024
dizziness   -0.008789153  0.01260153
### Visualization
shapiro.test(df$breathlessness)# NOT normally distributed

    Shapiro-Wilk normality test

data:  df$breathlessness
W = 0.80911, p-value = 2.605e-09
shapiro.test(df$dizziness)# NOT normally distributed

    Shapiro-Wilk normality test

data:  df$dizziness
W = 0.87343, p-value = 4.003e-07
shapiro.test(df$asthma)# NOT normally distributed

    Shapiro-Wilk normality test

data:  df$asthma
W = 0.51855, p-value = 1.662e-15
shapiro.test(df$unpleasantness)# NOT normally distributed

    Shapiro-Wilk normality test

data:  df$unpleasantness
W = 0.95164, p-value = 0.002433
discomfort1 <- jitteredViolin(df, df$unpleasantness, "A) Unpleasantness") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
  ylim(0, 120)
  
discomfort2 <- jitteredViolin(df, df$dizziness, "B) Dizziness") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
  ylim(0, 120)

discomfort3 <- jitteredViolin(df, df$breathlessness, "C) Breathlessness") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
  ylim(0, 120)

discomfort4 <- jitteredViolin(df, df$asthma, "D) Asthma") +
  stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1, label.y = 110) +
  ylim(0, 120)


gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2, nrow = 2)

  pdf("3.2_Plots_discomfort_violin.pdf")
  gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2)
  dev.off()
quartz_off_screen 
                2 
ggsave("Figure6.tiff", 
       plot =   gridExtra::grid.arrange(discomfort1, discomfort2, discomfort3, discomfort4, ncol = 2), 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

3.3) Response-Time

Note: for analysis on Response Time P029 had to be excluded as this patient was not able to press the response button themselves, but rather gave a sign for the experimenter to click with a mouseclick; thus the response time do not mirror properly the response time

# exclude particioant for all analysis on Resposne time and create new df
df_RT <- df[df$pcode != "29", ]

df_RT_HC <- df_RT[df_RT$group == "0", ]
df_RT_FND <- df_RT[df_RT$group == "1", ]

# mean age
df_RT %>% 
  group_by(group) %>% 
  get_summary_stats(average_RT, type = "mean_sd") 
# A tibble: 2 × 5
  group variable       n  mean    sd
  <fct> <fct>      <dbl> <dbl> <dbl>
1 0     average_RT    48  1.14 0.392
2 1     average_RT    42  1.07 0.414
median(df$average_RT[df$group=="1"])
[1] 0.98547
median(df$average_RT[df$group=="0"])
[1] 1.109339
IQR(df$average_RT[df$group=="1"])
[1] 0.5055542
IQR(df$average_RT[df$group=="0"])
[1] 0.5107342
res<-wilcox.test(df$average_RT[df$group=="1"], df$average_RT[df$group=="0"])
res # no sig difference between groups 

    Wilcoxon rank sum exact test

data:  df$average_RT[df$group == "1"] and df$average_RT[df$group == "0"]
W = 912, p-value = 0.3438
alternative hypothesis: true location shift is not equal to 0
## first we test whether response time correlates with the sensitivity score (as in, whether or not the sensitivity score is mirroring the impulsiveness or the jumping to conclusion)

cor.test(df_RT$average_RT, df_RT$rrst_sensitivity) # r = -0.26, p = 0.013

    Pearson's product-moment correlation

data:  df_RT$average_RT and df_RT$rrst_sensitivity
t = -2.5445, df = 88, p-value = 0.01269
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.44476592 -0.05783141
sample estimates:
       cor 
-0.2617873 
cor.test(df_RT$average_RT, df_RT$rrst_slope) # also testing the slope to represent interoceptive precision; p = 0.458

    Pearson's product-moment correlation

data:  df_RT$average_RT and df_RT$rrst_slope
t = 0.74444, df = 88, p-value = 0.4586
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1301134  0.2815872
sample estimates:
       cor 
0.07910932 
# separately per group
cor.test(df_RT_FND$average_RT, df_RT_FND$rrst_sensitivity) # p = 0.060

    Pearson's product-moment correlation

data:  df_RT_FND$average_RT and df_RT_FND$rrst_sensitivity
t = -1.9321, df = 40, p-value = 0.06045
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.54748026  0.01291418
sample estimates:
       cor 
-0.2921639 
cor.test(df_RT_FND$average_RT, df_RT_FND$rrst_slope) # p = 0.3745

    Pearson's product-moment correlation

data:  df_RT_FND$average_RT and df_RT_FND$rrst_slope
t = 0.99143, df = 40, p-value = 0.3274
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1564263  0.4381747
sample estimates:
      cor 
0.1548679 
cor.test(df_RT_HC$average_RT, df_RT_HC$rrst_sensitivity) # p = 0.041

    Pearson's product-moment correlation

data:  df_RT_HC$average_RT and df_RT_HC$rrst_sensitivity
t = -2.1061, df = 46, p-value = 0.04068
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.53556671 -0.01356923
sample estimates:
       cor 
-0.2965605 
cor.test(df_RT_HC$average_RT, df_RT_HC$rrst_slope) # p = 0.885

    Pearson's product-moment correlation

data:  df_RT_HC$average_RT and df_RT_HC$rrst_slope
t = 0.14545, df = 46, p-value = 0.885
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2643045  0.3037247
sample estimates:
       cor 
0.02144032 
# RRST Sensitivity GROUP controlled for average RT
rrst_rt <- lm(formula = rrst_sensitivity ~  group + average_RT,
                 data = df_RT) 
summary(rrst_rt)  # both sig

Call:
lm(formula = rrst_sensitivity ~ group + average_RT, data = df_RT)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.2230 -0.9557 -0.0913  1.0190  3.4978 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.6692     0.4472  10.441  < 2e-16 ***
group1       -0.7023     0.2834  -2.478  0.01514 *  
average_RT   -1.0012     0.3536  -2.831  0.00576 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.336 on 87 degrees of freedom
Multiple R-squared:  0.1299,    Adjusted R-squared:  0.1099 
F-statistic: 6.497 on 2 and 87 DF,  p-value: 0.002345
confint(rrst_rt, level = 0.95) 
                2.5 %     97.5 %
(Intercept)  3.780383  5.5580952
group1      -1.265640 -0.1390279
average_RT  -1.704135 -0.2983063
# Viszalization
df_plot <- df_RT %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))

RRST_RT<- ggplot(df_plot, aes(x = average_RT, y = rrst_sensitivity, color = as.factor(group))) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +  
  labs(x = "Average Response Time in s", 
       y = "Respiratory Sensitivity (RRST)", 
       title = "B) Respiratory Interoception and Response Time",
       color = "Group") +  
    theme_classic()+ 
  scale_color_manual(values = c("FND" = "#E95248", "HC" = "#868686FF")) +
  theme(panel.grid = element_blank())  
RRST_RT

ggsave("Figure5b.tiff", 
       plot = RRST_RT, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

RTplot <- plotViolinGroups_median(df_RT,df_RT$average_RT,"Average Response Time")+ stat_compare_means(method = "wilcox.test", paired = FALSE, label.x = 1.3)
RTplot

ggsave("Figure5a.tiff", 
       plot = RTplot, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

######
# Convert to long format
Data.RT <- select(df_RT, c("group", "incorrect_RT", "correct_RT", "pcode")) %>%
  mutate(group = factor(group, levels = c("0", "1"), labels = c("HC", "FND")))

df_long <- Data.RT %>%
  pivot_longer(cols = c(correct_RT, incorrect_RT),
               names_to = "accuracy",
               values_to = "RT") %>%
  mutate(
    accuracy = factor(accuracy, levels = c("correct_RT", "incorrect_RT"),
                      labels = c("Correct", "Incorrect")),
    group = factor(group, levels = c("FND", "HC")),
    group_accuracy = paste(group, accuracy, sep = "_")
  )

# Define color palette
group_colors <- c("FND" = "#E95248", "HC" = "#868686FF")
fill_colors <- c("FND_Correct" = "#f4cccc", 
                 "FND_Incorrect" = "#E95248",
                 "HC_Correct" = "#bcbcbc", 
                 "HC_Incorrect" = "#868686FF")

# Within-group comparisons (paired)
within_group_tests <- df_long %>%
  group_by(group) %>%
  wilcox_test(RT ~ accuracy, paired = TRUE) %>%
  add_significance() %>%
  mutate(group_accuracy = paste(group, "Correct_vs_Incorrect", sep = "_")) %>%
  add_xy_position(x = "group_accuracy", dodge = 0.6)

# Between-group comparisons (unpaired)
between_group_tests <- df_long %>%
  group_by(accuracy) %>%
  wilcox_test(RT ~ group, paired = FALSE) %>%
  add_significance() %>%
  mutate(group_accuracy = paste("FND", accuracy, sep = "_")) %>%  # use FND_* for x-position
  add_xy_position(x = "group_accuracy", dodge = 0.6)

violin_plots_RTs <- ggviolin(df_long, x = "group_accuracy", y = "RT",
         fill = "group_accuracy", palette = fill_colors,
         add = "jitter", width = 1, trim = FALSE) +
  labs(title = "A) Reaction Times by Group and Accuracy",
       x = NULL, y = "Reaction Time (RT)") +
  theme_classic() +
  theme(legend.position = "none")
violin_plots_RTs

ggsave("3.3_violin_plots_RTs.tiff", 
       plot = violin_plots_RTs,
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

ggsave("Figure4.tiff", 
       plot = gridExtra::grid.arrange(violin_plots_RTs, RRST_RT, ncol = 2), 
       dpi = 600, 
       width = 8, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

ggsave("Figure4a-b.svg", 
       plot = gridExtra::grid.arrange(violin_plots_RTs, RRST_RT, ncol = 2), 
       dpi = 600, 
       width = 8, 
       height = 4, 
       units = "in", 
       device = "svg")

4) Symptom Severity as dependent variables

# Compute the correlation matrix with variables of symptom severity to see whether they correlate
# Select the desired variables
severity_variables <- df_FND[, c("sfmdrs", "cgi", "sdq")]

# Compute the correlation matrix with rcorr
correlation_results <- Hmisc::rcorr(as.matrix(severity_variables))

# Extract correlation coefficients
cor_matrix <- correlation_results$r

# Extract p-values
p_matrix <- correlation_results$P

# Print results
cor_matrix
           sfmdrs       cgi        sdq
sfmdrs 1.00000000 0.5912098 0.05993439
cgi    0.59120983 1.0000000 0.30738677
sdq    0.05993439 0.3073868 1.00000000
p_matrix
             sfmdrs          cgi        sdq
sfmdrs           NA 2.986506e-05 0.70262560
cgi    2.986506e-05           NA 0.04495112
sdq    7.026256e-01 4.495112e-02         NA
# SIG correlations of different symptom severity scores
### CGI & SFMDRS corr: r= 0.5841928 , p=3.136863e-05
### CGI & SDQ corr: r = 0.31 , p = 0.045

4.1 ) Correlations with SDQ

# Exclude participants with sfmdrs = 0
filtered_df <- df_FND[df_FND$sfmdrs != 0, ]
# Count how many participants were excluded
excluded_count <- sum(df_FND$sfmdrs == 0)
excluded_count 
[1] 10
# SIMPLE CORRELATION OF RRST sensitivity and symptom severity scores - including SDQ
cor1 <- cor.test(df_FND$rrst_sensitivity, df_FND$sdq, method = "pearson") # r = -0.3828185, p-value = 0.01129
cor2 <- cor.test(filtered_df$rrst_sensitivity, filtered_df$sfmdrs, method = "pearson") 
cor3 <- cor.test(df_FND$rrst_sensitivity, df_FND$cgi, method = "pearson") 
p_values <- c(cor1$p.value, cor2$p.value, cor3$p.value) 
p_adjusted <- p.adjust(p_values, method = "bonferroni") 
results_rrst <- data.frame(
  Comparison = c("FND: rrst_sensitivity vs sdq", "FND: rrst_sensitivity vs sfmdrs", "FND: rrst_sensitivity vs cgi"),
  Original_P = p_values,
  Adjusted_P = p_adjusted)
print(results_rrst)
                       Comparison Original_P Adjusted_P
1    FND: rrst_sensitivity vs sdq 0.01128738 0.03386214
2 FND: rrst_sensitivity vs sfmdrs 0.36151691 1.00000000
3    FND: rrst_sensitivity vs cgi 0.30365832 0.91097497
# Save results as a text file
write.table(results_rrst, "3.1_corr_RRST_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)

cor.test(df$metascore, df$sdq, method = "pearson")  

    Pearson's product-moment correlation

data:  df$metascore and df$sdq
t = -2.6279, df = 89, p-value = 0.01012
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.44944344 -0.06604063
sample estimates:
       cor 
-0.2683359 
# p-value = 0.01012, cor = -0.2683359 

#  metascore
cor1 <- cor.test(df_FND$metascore, df_FND$sdq, method = "pearson")  # r = -0.36099, p-value = 0.0174
cor2 <- cor.test(filtered_df$metascore, filtered_df$sfmdrs, method = "pearson") 
cor3 <- cor.test(df_FND$metascore, df_FND$cgi, method = "pearson") 
p_values <- c(cor1$p.value, cor2$p.value, cor3$p.value) 
p_adjusted <- p.adjust(p_values, method = "bonferroni")  
results_meta <- data.frame(
  Comparison = c("FND: metascore vs sdq", "FND: metascore vs sfmdrs", "FND: metascore vs cgi"),
  Original_P = p_values,
  Adjusted_P = p_adjusted)
print(results_meta)
                Comparison Original_P Adjusted_P
1    FND: metascore vs sdq 0.01739608 0.05218824
2 FND: metascore vs sfmdrs 0.62780713 1.00000000
3    FND: metascore vs cgi 0.31527501 0.94582503
  # Save results as a text file
  write.table(results_rrst, "3.1_corr_meta_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)


#### CLINICAL VARIABLES and INTEROCEPTIVE SENSIBILITY
# separate per group
cor.test(df_FND$maia_total, df_FND$sdq, method = "pearson")  

    Pearson's product-moment correlation

data:  df_FND$maia_total and df_FND$sdq
t = -1.0818, df = 41, p-value = 0.2856
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4446866  0.1407959
sample estimates:
       cor 
-0.1665927 
cor.test(filtered_df$maia_total, filtered_df$sfmdrs, method = "pearson")  

    Pearson's product-moment correlation

data:  filtered_df$maia_total and filtered_df$sfmdrs
t = 1.1102, df = 31, p-value = 0.2754
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1583914  0.5049599
sample estimates:
      cor 
0.1955512 
cor.test(df_FND$maia_total, df_FND$cgi, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$maia_total and df_FND$cgi
t = 0.83381, df = 41, p-value = 0.4092
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1781232  0.4134383
sample estimates:
      cor 
0.1291288 
cor.test(df_FND$ias, df_FND$sdq, method = "pearson")  # t = -2.32, df = 41, p-value = 0.02539

    Pearson's product-moment correlation

data:  df_FND$ias and df_FND$sdq
t = -2.32, df = 41, p-value = 0.02539
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.58150466 -0.04490785
sample estimates:
      cor 
-0.340657 
cor.test(filtered_df$ias, filtered_df$sfmdrs, method = "pearson")  

    Pearson's product-moment correlation

data:  filtered_df$ias and filtered_df$sfmdrs
t = 0.6599, df = 31, p-value = 0.5142
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.2351120  0.4431019
sample estimates:
      cor 
0.1176971 
cor.test(df_FND$ias, df_FND$cgi, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$ias and df_FND$cgi
t = -1.602, df = 41, p-value = 0.1168
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.50615656  0.06216586
sample estimates:
       cor 
-0.2427097 
# Collect all p-values from test of interoception varibales and clincial variables
p_values <- c(
  cor.test(df_FND$maia_total, df_FND$sdq, method = "pearson")$p.value,
  cor.test(filtered_df$maia_total, filtered_df$sfmdrs, method = "pearson")$p.value,
  cor.test(df_FND$maia_total, df_FND$cgi, method = "pearson")$p.value,
  
  cor.test(df_FND$ias, df_FND$sdq, method = "pearson")$p.value,
  cor.test(filtered_df$ias, filtered_df$sfmdrs, method = "pearson")$p.value,
  cor.test(df_FND$ias, df_FND$cgi, method = "pearson")$p.value
)

# Apply Bonferroni correction
adjusted_p_values <- p.adjust(p_values, method = "bonferroni")

# Combine with test names for clarity
test_names <- c(
  "MAIA vs sdq (FND)",
  "MAIA vs sfmdrs (FND)",
  "MAIAvs cgi (FND)",
  
  "IAS vs sdq (FND)",
  "IAS vs sfmdrs (FND)",
  "IAS vs cgi (FND)")

# Create a data frame for easier interpretation
results <- data.frame(
  Test = test_names,
  P_Value = p_values,
  Adjusted_P_Value = adjusted_p_values)

# View results
print(results)
                  Test   P_Value Adjusted_P_Value
1    MAIA vs sdq (FND) 0.2856496        1.0000000
2 MAIA vs sfmdrs (FND) 0.2754433        1.0000000
3     MAIAvs cgi (FND) 0.4092203        1.0000000
4     IAS vs sdq (FND) 0.0253949        0.1523694
5  IAS vs sfmdrs (FND) 0.5141939        1.0000000
6     IAS vs cgi (FND) 0.1168329        0.7009974
  # Save results as a text file
  write.table(results_rrst, "3.1_corr_intero_symptoms.txt", sep = "\t", row.names = FALSE, quote = FALSE)

4.2 ) Identifying Covariates for Regression Analysis

We identified SDQ to be correlated with RRST sensitivity and metacognition in the FND population only.

Next, we need to identify which covariates will be included in our linear regression (separate per group) to test whether these correlation survive also the control for these covariates. To do so, we choose those variables that are both theoretically but also statistically correlated with both the predictor (RRST sensitivity and metacognition in our case) and the outcome variable (SDQ in our case).

# does age need to be a covariate in the model ?
cor.test(df_FND$age, df_FND$rrst_sensitivity, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$age and df_FND$rrst_sensitivity
t = -0.99609, df = 41, p-value = 0.325
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4340205  0.1537269
sample estimates:
      cor 
-0.153714 
cor.test(df_FND$age, df_FND$metascore, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$age and df_FND$metascore
t = 1.0698, df = 41, p-value = 0.291
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1426087  0.4432013
sample estimates:
      cor 
0.1647936 
cor.test(df_FND$age, df_FND$sdq, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$age and df_FND$sdq
t = 1.6267, df = 41, p-value = 0.1115
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.05844549  0.50892864
sample estimates:
      cor 
0.2462205 
# does bmi need to be a covariate in the model?
cor.test(df_FND$bmi, df_FND$rrst_sensitivity, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$bmi and df_FND$rrst_sensitivity
t = -0.5235, df = 41, p-value = 0.6034
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3727075  0.2243494
sample estimates:
        cor 
-0.08148511 
cor.test(df_FND$bmi, df_FND$metascore, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$bmi and df_FND$metascore
t = -0.58658, df = 41, p-value = 0.5607
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3811282  0.2150079
sample estimates:
        cor 
-0.09122709 
cor.test(df_FND$bmi, df_FND$sdq, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$bmi and df_FND$sdq
t = -0.97486, df = 41, p-value = 0.3353
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4313571  0.1569247
sample estimates:
       cor 
-0.1505131 
# does anx_dep_sum need to be a covariate in the model?
cor.test(df_FND$anx_dep_SUM, df_FND$rrst_sensitivity, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$anx_dep_SUM and df_FND$rrst_sensitivity
t = -1.0783, df = 41, p-value = 0.2872
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4442492  0.1413301
sample estimates:
       cor 
-0.1660627 
cor.test(df_FND$anx_dep_SUM, df_FND$metascore, method = "pearson") 

    Pearson's product-moment correlation

data:  df_FND$anx_dep_SUM and df_FND$metascore
t = -0.64011, df = 41, p-value = 0.5257
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.3882177  0.2070577
sample estimates:
       cor 
-0.0994722 
cor.test(df_FND$anx_dep_SUM, df_FND$sdq, method = "pearson") # trend with r= 0.294215, p-value = 0.05548

    Pearson's product-moment correlation

data:  df_FND$anx_dep_SUM and df_FND$sdq
t = 1.9711, df = 41, p-value = 0.05548
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.006722943  0.546285935
sample estimates:
     cor 
0.294215 
# does medication need to be a covariate in the model?
res<-t.test(df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "1"], df_FND$rrst_sensitivity[df_FND$psychotropic_medication=="0"])
res # sig difference across groups with or without intake of psychotropic medicatioN: t = -2.1273, p-value = 0.04011

    Welch Two Sample t-test

data:  df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "1"] and df_FND$rrst_sensitivity[df_FND$psychotropic_medication == "0"]
t = -2.1273, df = 37.077, p-value = 0.04011
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.00230723 -0.04882873
sample estimates:
mean of x mean of y 
 2.309574  3.335142 
res<-t.test(df_FND$metascore[df_FND$psychotropic_medication == "1"], df_FND$metascore[df_FND$psychotropic_medication=="0"])
res # 

    Welch Two Sample t-test

data:  df_FND$metascore[df_FND$psychotropic_medication == "1"] and df_FND$metascore[df_FND$psychotropic_medication == "0"]
t = -1.4203, df = 39.126, p-value = 0.1634
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.07204124  0.01260121
sample estimates:
mean of x mean of y 
0.5533079 0.5830279 
res<-t.test(df_FND$sdq[df_FND$psychotropic_medication == "1"], df_FND$sdq[df_FND$psychotropic_medication=="0"])
res # sig difference across groups with or without intake of psychotropic medication: t = 2.6337, p-value = 0.0129

    Welch Two Sample t-test

data:  df_FND$sdq[df_FND$psychotropic_medication == "1"] and df_FND$sdq[df_FND$psychotropic_medication == "0"]
t = 2.6337, df = 32.012, p-value = 0.0129
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  2.584876 20.230913
sample estimates:
mean of x mean of y 
 46.15789  34.75000 
# differences of sex in our predictor and outcome variable?
res<-t.test(df_FND$rrst_sensitivity[df_FND$sex=="1"], df_FND$rrst_sensitivity[df_FND$sex=="2"])
res # 

    Welch Two Sample t-test

data:  df_FND$rrst_sensitivity[df_FND$sex == "1"] and df_FND$rrst_sensitivity[df_FND$sex == "2"]
t = -1.1197, df = 15.303, p-value = 0.2801
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.0795758  0.6455268
sample estimates:
mean of x mean of y 
 2.681884  3.398908 
res<-t.test(df_FND$metascore[df_FND$sex=="1"], df_FND$metascore[df_FND$sex=="2"])
res # 

    Welch Two Sample t-test

data:  df_FND$metascore[df_FND$sex == "1"] and df_FND$metascore[df_FND$sex == "2"]
t = 1.1436, df = 16.994, p-value = 0.2687
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.02470859  0.08319037
sample estimates:
mean of x mean of y 
0.5780561 0.5488152 
res<-t.test(df_FND$sdq[df_FND$sex=="1"], df_FND$sdq[df_FND$sex=="2"])
res # 

    Welch Two Sample t-test

data:  df_FND$sdq[df_FND$sex == "1"] and df_FND$sdq[df_FND$sex == "2"]
t = 0.33612, df = 17.95, p-value = 0.7407
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.402255 12.982900
sample estimates:
mean of x mean of y 
 40.29032  38.50000 

Seeing that the intake of psychotropic medication is both correlated to our outcome and predictor variables, we have to include this in our linear regression. We add it in an interaction term with group, as it is correlated with group, and then we run it on the full model (including both groups). For the other variables (sex, age, bmi) we do not have to include those, as they are not statistically correlated to both variables. For affective symptoms - there is a trend, hence we run a supplemantary analysis for including these as well.

4.3) Linear Regression

# RRST sensitivity
sdq_lm <- lm(formula = sdq ~  rrst_sensitivity*group + psychotropic_medication,
                 data = df) 
summary(sdq_lm)  

Call:
lm(formula = sdq ~ rrst_sensitivity * group + psychotropic_medication, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.741  -4.157  -1.555   3.758  38.900 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               22.4719     4.8244   4.658 1.18e-05 ***
rrst_sensitivity           0.3626     1.2835   0.282 0.778274    
group1                    21.9827     5.8059   3.786 0.000286 ***
psychotropic_medication1   7.2866     2.8532   2.554 0.012460 *  
rrst_sensitivity:group1   -3.0980     1.5711  -1.972 0.051916 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.575 on 84 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.4898,    Adjusted R-squared:  0.4655 
F-statistic: 20.16 on 4 and 84 DF,  p-value: 1.146e-11
confint(sdq_lm, level = 0.95)  
                             2.5 %      97.5 %
(Intercept)              12.878160 32.06567055
rrst_sensitivity         -2.189855  2.91497870
group1                   10.436956 33.52837270
psychotropic_medication1  1.612716 12.96052136
rrst_sensitivity:group1  -6.222352  0.02630844
  # Save as a text file
  lm_summary <- capture.output(summary(sdq_lm))
  lm_confint <- capture.output(confint(sdq_lm, level = 0.95))
  lm_results <- c("### Linear Model Summary ###", lm_summary, 
                "", "### Confidence Intervals ###", lm_confint)
  writeLines(lm_results, "3.3_lm_rrst_med.txt")


sdq_lm_affective <- lm(formula = sdq ~  rrst_sensitivity*group + psychotropic_medication + anx_dep_SUM,
                 data = df) 
summary(sdq_lm_affective)  # 

Call:
lm(formula = sdq ~ rrst_sensitivity * group + psychotropic_medication + 
    anx_dep_SUM, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.061  -4.445  -1.013   3.018  38.181 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              19.57599    5.06195   3.867 0.000218 ***
rrst_sensitivity         -0.03933    1.29067  -0.030 0.975765    
group1                   19.12577    5.97891   3.199 0.001954 ** 
psychotropic_medication1  5.04848    3.11010   1.623 0.108327    
anx_dep_SUM               0.10520    0.06154   1.710 0.091094 .  
rrst_sensitivity:group1  -2.67953    1.57261  -1.704 0.092146 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.468 on 83 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.5071,    Adjusted R-squared:  0.4775 
F-statistic: 17.08 on 5 and 83 DF,  p-value: 1.39e-11
confint(sdq_lm_affective, level = 0.95)  
                               2.5 %     97.5 %
(Intercept)               9.50798058 29.6440034
rrst_sensitivity         -2.60641766  2.5277619
group1                    7.23395740 31.0175789
psychotropic_medication1 -1.13737894 11.2343361
anx_dep_SUM              -0.01719769  0.2276056
rrst_sensitivity:group1  -5.80738085  0.4483295
  # Save as a text file
  lm_summary <- capture.output(summary(sdq_lm_affective))
  lm_confint <- capture.output(confint(sdq_lm_affective, level = 0.95))
  lm_results <- c("### Linear Model Summary ###", lm_summary, 
                "", "### Confidence Intervals ###", lm_confint)
  writeLines(lm_results, "3.3_lm_rrst_med_aff.txt")


# metascore
sdq_lm <- lm(formula = sdq ~  metascore*group + psychotropic_medication,
                 data = df) 
summary(sdq_lm)  # 

Call:
lm(formula = sdq ~ metascore * group + psychotropic_medication, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.427  -3.991  -1.279   3.525  38.022 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
(Intercept)                18.253     21.722   0.840   0.4031   
metascore                   9.522     37.609   0.253   0.8007   
group1                     54.342     25.040   2.170   0.0328 * 
psychotropic_medication1    7.997      2.760   2.897   0.0048 **
metascore:group1          -73.283     43.279  -1.693   0.0941 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.56 on 84 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.4914,    Adjusted R-squared:  0.4672 
F-statistic: 20.29 on 4 and 84 DF,  p-value: 1.005e-11
confint(sdq_lm, level = 0.95)  
                               2.5 %    97.5 %
(Intercept)               -24.944466  61.44956
metascore                 -65.266950  84.31116
group1                      4.546173 104.13700
psychotropic_medication1    2.508047  13.48622
metascore:group1         -159.348639  12.78268
  # Save as a text file
  lm_summary <- capture.output(summary(sdq_lm))
  lm_confint <- capture.output(confint(sdq_lm, level = 0.95))
  lm_results <- c("### Linear Model Summary ###", lm_summary, 
                "", "### Confidence Intervals ###", lm_confint)
  writeLines(lm_results, "3.3_lm_meta_med.txt")

sdq_lm_affective <- lm(formula = sdq ~  metascore*group + psychotropic_medication + anx_dep_SUM,
                 data = df) 
summary(sdq_lm_affective)  # 

Call:
lm(formula = sdq ~ metascore * group + psychotropic_medication + 
    anx_dep_SUM, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.280  -4.199  -1.335   3.362  34.590 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)  
(Intercept)               15.41127   21.50568   0.717   0.4756  
metascore                  6.75580   37.16375   0.182   0.8562  
group1                    51.34935   24.77943   2.072   0.0413 *
psychotropic_medication1   5.78845    2.99396   1.933   0.0566 .
anx_dep_SUM                0.10755    0.06037   1.782   0.0785 .
metascore:group1         -70.57624   42.75691  -1.651   0.1026  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.439 on 83 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.5102,    Adjusted R-squared:  0.4806 
F-statistic: 17.29 on 5 and 83 DF,  p-value: 1.088e-11
confint(sdq_lm_affective, level = 0.95)  
                                 2.5 %      97.5 %
(Intercept)               -27.36266743  58.1852148
metascore                 -67.16139890  80.6729973
group1                      2.06405554 100.6346377
psychotropic_medication1   -0.16641925  11.7433184
anx_dep_SUM                -0.01252055   0.2276273
metascore:group1         -155.61801619  14.4655269
  # Save as a text file
  lm_summary <- capture.output(summary(sdq_lm_affective))
  lm_confint <- capture.output(confint(sdq_lm_affective, level = 0.95))
  lm_results <- c("### Linear Model Summary ###", lm_summary, 
                "", "### Confidence Intervals ###", lm_confint)
  writeLines(lm_results, "3.3_lm_meta_med_aff.txt")

4.4) SUPPLEMENTARY Regression - separate per group and controlling for ANX-DEP sum scores

# RRST sensitivity FND
sdq_lm_FND <- lm(formula = sdq ~  rrst_sensitivity + psychotropic_medication,
                 data = df_FND) 
summary(sdq_lm_FND)  

Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication, 
    data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.141  -8.096  -2.786   6.417  37.778 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)                43.391      5.165   8.401 2.28e-10 ***
rrst_sensitivity           -2.591      1.323  -1.959   0.0571 .  
psychotropic_medication1    8.751      4.263   2.053   0.0467 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.16 on 40 degrees of freedom
Multiple R-squared:  0.2279,    Adjusted R-squared:  0.1893 
F-statistic: 5.903 on 2 and 40 DF,  p-value: 0.005668
confint(sdq_lm_FND, level = 0.95)  
                              2.5 %      97.5 %
(Intercept)              32.9527048 53.82919851
rrst_sensitivity         -5.2640573  0.08229745
psychotropic_medication1  0.1358901 17.36565231
sdq_lm_affective_FND <- lm(formula = sdq ~  rrst_sensitivity + psychotropic_medication + anx_dep_SUM,
                 data = df_FND) 
summary(sdq_lm_affective_FND)  # 

Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication + 
    anx_dep_SUM, data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.312  -8.076  -2.783   5.705  37.430 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              39.34986    7.81643   5.034 1.12e-05 ***
rrst_sensitivity         -2.60166    1.33144  -1.954   0.0579 .  
psychotropic_medication1  6.88075    5.06997   1.357   0.1825    
anx_dep_SUM               0.07699    0.11120   0.692   0.4928    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.25 on 39 degrees of freedom
Multiple R-squared:  0.2373,    Adjusted R-squared:  0.1786 
F-statistic: 4.044 on 3 and 39 DF,  p-value: 0.0135
confint(sdq_lm_affective_FND, level = 0.95)  
                              2.5 %      97.5 %
(Intercept)              23.5396403 55.16007691
rrst_sensitivity         -5.2947400  0.09142542
psychotropic_medication1 -3.3742320 17.13572965
anx_dep_SUM              -0.1479358  0.30190731
# RRST sensitivity HC
sdq_lm_HC <- lm(formula = sdq ~  rrst_sensitivity + psychotropic_medication,
                 data = df_HC) 
summary(sdq_lm_HC)  

Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication, 
    data = df_HC)

Residuals:
   Min     1Q Median     3Q    Max 
-4.328 -3.163 -1.072  1.924 11.822 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              24.45286    2.05771  11.884 3.57e-15 ***
rrst_sensitivity         -0.09285    0.54018  -0.172    0.864    
psychotropic_medication1 -0.76881    2.93689  -0.262    0.795    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.865 on 43 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.001805,  Adjusted R-squared:  -0.04462 
F-statistic: 0.03888 on 2 and 43 DF,  p-value: 0.9619
confint(sdq_lm_HC, level = 0.95)  
                             2.5 %     97.5 %
(Intercept)              20.303100 28.6026252
rrst_sensitivity         -1.182222  0.9965168
psychotropic_medication1 -6.691613  5.1539933
sdq_lm_affective_HC <- lm(formula = sdq ~  rrst_sensitivity + psychotropic_medication + anx_dep_SUM,
                 data = df_HC) 
summary(sdq_lm_affective_HC)  # 

Call:
lm(formula = sdq ~ rrst_sensitivity + psychotropic_medication + 
    anx_dep_SUM, data = df_HC)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.1591 -2.8961 -0.4383  1.9792  8.4374 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              20.46811    2.17139   9.426 6.35e-12 ***
rrst_sensitivity         -0.45618    0.49392  -0.924   0.3610    
psychotropic_medication1 -1.35880    2.62913  -0.517   0.6080    
anx_dep_SUM               0.12607    0.03657   3.447   0.0013 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.453 on 42 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.222, Adjusted R-squared:  0.1664 
F-statistic: 3.994 on 3 and 42 DF,  p-value: 0.0137
confint(sdq_lm_affective_HC, level = 0.95)  
                               2.5 %     97.5 %
(Intercept)              16.08607859 24.8501447
rrst_sensitivity         -1.45294757  0.5405941
psychotropic_medication1 -6.66459477  3.9470002
anx_dep_SUM               0.05226925  0.1998698
# metascore FND
sdq_lm_FND <- lm(formula = sdq ~  metascore + psychotropic_medication,
                 data = df_FND) 
summary(sdq_lm_FND)  # 

Call:
lm(formula = sdq ~ metascore + psychotropic_medication, data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.720  -7.438  -3.396   8.541  37.124 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)                70.469     17.675   3.987 0.000277 ***
metascore                 -61.264     29.966  -2.044 0.047528 *  
psychotropic_medication1    9.587      4.123   2.326 0.025200 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.11 on 40 degrees of freedom
Multiple R-squared:  0.2339,    Adjusted R-squared:  0.1956 
F-statistic: 6.106 on 2 and 40 DF,  p-value: 0.004851
confint(sdq_lm_FND, level = 0.95)  
                               2.5 %     97.5 %
(Intercept)                34.747037 106.190457
metascore                -121.827371  -0.701054
psychotropic_medication1    1.255103  17.919140
sdq_lm_affective_FND <- lm(formula = sdq ~  metascore + psychotropic_medication + anx_dep_SUM,
                 data = df_FND) 
summary(sdq_lm_affective_FND)  # 

Call:
lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM, 
    data = df_FND)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.069  -7.355  -2.813   7.833  34.746 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)               66.54227   18.59586   3.578 0.000944 ***
metascore                -61.78579   30.15574  -2.049 0.047244 *  
psychotropic_medication1   7.64265    4.94639   1.545 0.130400    
anx_dep_SUM                0.07988    0.11073   0.721 0.474958    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 13.19 on 39 degrees of freedom
Multiple R-squared:  0.244, Adjusted R-squared:  0.1858 
F-statistic: 4.195 on 3 and 39 DF,  p-value: 0.0115
confint(sdq_lm_affective_FND, level = 0.95)  
                                2.5 %      97.5 %
(Intercept)                28.9286047 104.1559436
metascore                -122.7815257  -0.7900468
psychotropic_medication1   -2.3623702  17.6476663
anx_dep_SUM                -0.1440951   0.3038646
# metascore HC

sdq_lm_HC <- lm(formula = sdq ~  metascore + psychotropic_medication,
                 data = df_HC) 
summary(sdq_lm_HC)  # 

Call:
lm(formula = sdq ~ metascore + psychotropic_medication, data = df_HC)

Residuals:
   Min     1Q Median     3Q    Max 
-4.280 -3.087 -1.230  2.018 11.924 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)  
(Intercept)               21.7955     8.8443   2.464   0.0178 *
metascore                  4.0164    15.2906   0.263   0.7941  
psychotropic_medication1  -0.5258     2.8134  -0.187   0.8526  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.864 on 43 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.002719,  Adjusted R-squared:  -0.04367 
F-statistic: 0.05863 on 2 and 43 DF,  p-value: 0.9431
confint(sdq_lm_HC, level = 0.95)  
                              2.5 %    97.5 %
(Intercept)                3.959307 39.631774
metascore                -26.820060 34.852853
psychotropic_medication1  -6.199508  5.147904
sdq_lm_affective_HC <- lm(formula = sdq ~  metascore + psychotropic_medication + anx_dep_SUM,
                 data = df_HC) 
summary(sdq_lm_affective_HC)  # 

Call:
lm(formula = sdq ~ metascore + psychotropic_medication + anx_dep_SUM, 
    data = df_HC)

Residuals:
   Min     1Q Median     3Q    Max 
-5.935 -2.825 -0.655  2.306  8.960 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)   
(Intercept)              17.65519    8.07978   2.185  0.03451 * 
metascore                 2.52815   13.80551   0.183  0.85558   
psychotropic_medication1 -0.54217    2.53876  -0.214  0.83193   
anx_dep_SUM               0.11865    0.03609   3.287  0.00205 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.487 on 42 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.2068,    Adjusted R-squared:  0.1501 
F-statistic:  3.65 on 3 and 42 DF,  p-value: 0.01995
confint(sdq_lm_affective_HC, level = 0.95)  
                               2.5 %     97.5 %
(Intercept)                1.3495342 33.9608409
metascore                -25.3324965 30.3887882
psychotropic_medication1  -5.6655976  4.5812588
anx_dep_SUM                0.0458067  0.1914852

5) Modeling Accuracy Interaction

R-script provided from Reviewer 3 and adapted to our data; thus acknolwedgment to them

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0     ✔ tibble  3.2.1
✔ purrr   1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ rstatix::filter()  masks dplyr::filter(), stats::filter()
✖ dplyr::lag()       masks stats::lag()
✖ car::recode()      masks dplyr::recode()
✖ purrr::some()      masks car::some()
✖ Hmisc::src()       masks dplyr::src()
✖ Hmisc::summarize() masks dplyr::summarize()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(glmmTMB)
Warning in check_dep_version(dep_pkg = "TMB"): package version mismatch: 
glmmTMB was built with TMB package version 1.9.17
Current TMB package version is 1.9.10
Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
library(ggeffects)
library(patchwork)
library(sjPlot)
Learn more about sjPlot with 'browseVignettes("sjPlot")'.
library(ggplot2)
library(ggthemes)

# Example RRST analysis script: modeling confidence data
# Assumes 'rrst_data' contains:
# - confidence: trial-level confidence ratings (scaled 0-1)
# - stimLevel: stimulus intensity (continuous or ordinal)
# - ResponseAccuracy: binary (0 = incorrect, 1 = correct)
# - group: factor (e.g., "FND", "Control")
# - SDQ_score: continuous clinical variable (e.g., dissociation)    
# - pcode: participant ID   

rrst_confidence_analysis <- readxl::read_xlsx( "/Users/nataschastoffel/Documents/GitHub/interoception_NS/data/processed/05-2025/rrst_confidence_analysis.xlsx", na = c("NaN", "NA", ""))
# data is in wide format, so we have to change to long

library(tidyverse)

rrst_data <- rrst_confidence_analysis %>%
  pivot_longer(
    cols = matches("^(accuracy|stim|conf)_trial\\d+$"),
    names_to = c(".value", "trial"),
    names_pattern = "(accuracy|stim|conf)_trial(\\d+)"
  ) %>%
  rename(
    ResponseAccuracy = accuracy,
    stimLevel = stim,
    confidence = conf
  ) %>%
  mutate(
    trial = as.integer(trial),
    ResponseAccuracy = as.integer(ResponseAccuracy),  # binary
    stimLevel = as.numeric(stimLevel),                # continuous
    confidence = as.numeric(confidence),              # continuous
    group = as.factor(group),
    sdq = as.numeric(sdq),
    pcode = as.character(pcode) )


# Note: Consider including (1 | stimLevel) as a random effect to control for stimulus-level variance if needed

rrst_data <- rrst_data %>%
  mutate(
    ResponseAccuracy = factor(ResponseAccuracy, levels = c(1, 0), labels = c("Correct", "Incorrect")),
    group = factor(group),  # Make sure group is also a factor if needed
    pcode = as.character(pcode)
  )

5.1) Interaction with GROUP

# Model: Confidence ~ Stimulus × Accuracy × Group + (1|pcode)
conf_fit <- glmmTMB(confidence ~ stimLevel * ResponseAccuracy * group + (1 | pcode) + (1 | stimLevel),
                   data = rrst_data,
                   family = ordbeta(),
                   start = list(psi = c(0, 1)))

summary(conf_fit)
 Family: ordbeta  ( logit )
Formula:          
confidence ~ stimLevel * ResponseAccuracy * group + (1 | pcode) +  
    (1 | stimLevel)
Data: rrst_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   3970.6    4053.9   -1972.3    3944.6      4470 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 pcode     (Intercept) 0.37135  0.6094  
 stimLevel (Intercept) 0.04017  0.2004  
Number of obs: 4483, groups:  pcode, 91; stimLevel, 17

Dispersion parameter for ordbeta family (): 3.19 

Conditional model:
                                            Estimate Std. Error z value
(Intercept)                                  0.22137    0.29777   0.743
stimLevel                                    0.03737    0.02184   1.711
ResponseAccuracyIncorrect                    0.10508    0.28360   0.370
groupHC                                     -0.42289    0.30412  -1.391
stimLevel:ResponseAccuracyIncorrect         -0.04768    0.02088  -2.283
stimLevel:groupHC                            0.03203    0.02027   1.580
ResponseAccuracyIncorrect:groupHC            0.26384    0.36591   0.721
stimLevel:ResponseAccuracyIncorrect:groupHC -0.01864    0.02769  -0.673
                                            Pr(>|z|)  
(Intercept)                                   0.4572  
stimLevel                                     0.0871 .
ResponseAccuracyIncorrect                     0.7110  
groupHC                                       0.1644  
stimLevel:ResponseAccuracyIncorrect           0.0224 *
stimLevel:groupHC                             0.1141  
ResponseAccuracyIncorrect:groupHC             0.4709  
stimLevel:ResponseAccuracyIncorrect:groupHC   0.5008  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model summary table
tab_model(conf_fit,
         pred.labels = c("Intercept",
                         "Stimulus Intensity",
                         "Response Accuracy",
                         "Group (FND vs Control)",
                         "Stim × Accuracy",
                         "Stim × Group",
                         "Accuracy × Group",
                         "Stim × Accuracy × Group"),
         dv.labels = "Confidence",
         file = "RRST_confidence_model.doc")
  Confidence
Predictors Estimates CI p
Intercept 1.25 0.70 – 2.24 0.457
Stimulus Intensity 1.04 0.99 – 1.08 0.087
Response Accuracy 1.11 0.64 – 1.94 0.711
Group (FND vs Control) 0.66 0.36 – 1.19 0.164
Stim × Accuracy 0.95 0.92 – 0.99 0.022
Stim × Group 1.03 0.99 – 1.07 0.114
Accuracy × Group 1.30 0.64 – 2.67 0.471
Stim × Accuracy × Group 0.98 0.93 – 1.04 0.501
Random Effects
σ2 0.14
τ00 pcode 0.37
τ00 stimLevel 0.04
ICC 0.75
N pcode 91
N stimLevel 17
Observations 4483
Marginal R2 / Conditional R2 0.106 / 0.775
# Define palette
color_palette <- c("Correct" = "#1b9e77", "Incorrect" = "#d95f02")

# Plot 1: Group × Accuracy Interaction
plot_group_accuracy <- plot_model(conf_fit, type = "pred", terms = c("group", "ResponseAccuracy")) +
 scale_color_manual(values = color_palette) +
 labs(title = "Group and Accuracy Interaction",
      x = "Group",
      y = "Predicted Confidence") +
 theme_minimal(base_size = 14) +
 theme(legend.position = "top",
       plot.title = element_text(face = "bold", hjust = 0.5),
       axis.title = element_text(size = 14),
       axis.text = element_text(size = 12))

# Plot 2: Stimulus × Accuracy × Group Interaction
plot_stim_acc_group <- plot_model(conf_fit, type = "pred", terms = c("stimLevel", "ResponseAccuracy", "group")) +
  scale_color_manual(values = color_palette) +
  scale_fill_manual(values = color_palette) + 
  labs(title = "Stimulus Intensity × Accuracy × Group Interaction",
       x = "Stimulus Intensity",
       y = "Predicted Confidence") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top",
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12))

# Combine plots
combined_plot <- plot_group_accuracy / plot_stim_acc_group +
 plot_layout(guides = 'collect') & theme(legend.position = "bottom")

print(combined_plot)

ggsave("S2_Accuracy_interactions_GROUP.tiff", 
       plot = combined_plot, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")

5.2) Interaction with SDQ

# Model: Confidence ~ Stimulus × Accuracy × Group + (1|pcode)
conf_fit <- glmmTMB(confidence ~ stimLevel * ResponseAccuracy * sdq + (1 | pcode) + (1 | stimLevel),
                   data = rrst_data,
                   family = ordbeta(),
                   start = list(psi = c(0, 1)))
summary(conf_fit)
 Family: ordbeta  ( logit )
Formula:          
confidence ~ stimLevel * ResponseAccuracy * sdq + (1 | pcode) +  
    (1 | stimLevel)
Data: rrst_data

      AIC       BIC    logLik -2*log(L)  df.resid 
   3930.6    4013.9   -1952.3    3904.6      4470 

Random effects:

Conditional model:
 Groups    Name        Variance Std.Dev.
 pcode     (Intercept) 0.36705  0.6058  
 stimLevel (Intercept) 0.04832  0.2198  
Number of obs: 4483, groups:  pcode, 91; stimLevel, 17

Dispersion parameter for ordbeta family (): 3.22 

Conditional model:
                                          Estimate Std. Error z value Pr(>|z|)
(Intercept)                             -1.0526862  0.4103804  -2.565  0.01031
stimLevel                                0.1552625  0.0288582   5.380 7.44e-08
ResponseAccuracyIncorrect                1.1060807  0.4739671   2.334  0.01961
sdq                                      0.0331362  0.0115360   2.872  0.00407
stimLevel:ResponseAccuracyIncorrect     -0.1499379  0.0369706  -4.056 5.00e-05
stimLevel:sdq                           -0.0031433  0.0007546  -4.165 3.11e-05
ResponseAccuracyIncorrect:sdq           -0.0232996  0.0142408  -1.636  0.10181
stimLevel:ResponseAccuracyIncorrect:sdq  0.0025795  0.0010877   2.372  0.01772
                                           
(Intercept)                             *  
stimLevel                               ***
ResponseAccuracyIncorrect               *  
sdq                                     ** 
stimLevel:ResponseAccuracyIncorrect     ***
stimLevel:sdq                           ***
ResponseAccuracyIncorrect:sdq              
stimLevel:ResponseAccuracyIncorrect:sdq *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Model summary table
tab_model(conf_fit,
         pred.labels = c("Intercept",
                         "Stimulus Intensity",
                         "Response Accuracy",
                         "SDQ",
                         "Stim × Accuracy",
                         "Stim × SDQ",
                         "Accuracy × SDQ",
                         "Stim × Accuracy × SDQ"),
         dv.labels = "Confidence",
         file = "RRST_confidence_model.doc")
  Confidence
Predictors Estimates CI p
Intercept 0.35 0.16 – 0.78 0.010
Stimulus Intensity 1.17 1.10 – 1.24 <0.001
Response Accuracy 3.02 1.19 – 7.65 0.020
SDQ 1.03 1.01 – 1.06 0.004
Stim × Accuracy 0.86 0.80 – 0.93 <0.001
Stim × SDQ 1.00 1.00 – 1.00 <0.001
Accuracy × SDQ 0.98 0.95 – 1.00 0.102
Stim × Accuracy × SDQ 1.00 1.00 – 1.00 0.018
Random Effects
σ2 0.14
τ00 pcode 0.37
τ00 stimLevel 0.05
ICC 0.75
N pcode 91
N stimLevel 17
Observations 4483
Marginal R2 / Conditional R2 0.129 / 0.784
# Define palette
color_palette <- c("Correct" = "#1b9e77", "Incorrect" = "#d95f02")

# Plot 1: SDQ × Accuracy Interaction
plot_sdq_accuracy <- plot_model(conf_fit, type = "pred", terms = c("sdq", "ResponseAccuracy")) +
 scale_color_manual(values = color_palette) +
 labs(title = "SDQ and Accuracy Interaction",
      x = "SDQ",
      y = "Predicted Confidence") +
 theme_minimal(base_size = 14) +
 theme(legend.position = "top",
       plot.title = element_text(face = "bold", hjust = 0.5),
       axis.title = element_text(size = 14),
       axis.text = element_text(size = 12))

# Plot 2: Stimulus × Accuracy × SDQ Interaction
plot_stim_acc_sdq <- plot_model(conf_fit, type = "pred", terms = c("stimLevel", "ResponseAccuracy", "sdq")) +
  scale_color_manual(values = color_palette) +
  scale_fill_manual(values = color_palette) + 
  labs(title = "Stimulus Intensity × Accuracy × SDQ Interaction",
       x = "Stimulus Intensity",
       y = "Predicted Confidence") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "top",
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 12))

# Combine plots
combined_plot <- plot_sdq_accuracy / plot_stim_acc_sdq +
 plot_layout(guides = 'collect') & theme(legend.position = "bottom")

print(combined_plot)

ggsave("S3_Accuracy_interactions_SDQ.tiff", 
       plot = combined_plot, 
       dpi = 600, 
       width = 6, 
       height = 4, 
       units = "in", 
       device = "tiff", 
       compression = "lzw")